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FOREWORD 


Mathematical  models,  methods  and  techniques  which  are  useful  and  appropriate  for  estimating 
the  accuracy  of  position  (in  M>me  coordinate  systems);  velocity,  and  acceleration  data  are  presented  in 
this  document.  The  development  and  use  of  the  techniques  discussed  have  evolved  through  the  years 
and  in  some  cases  out  of  work  not  related  to  missile  testing.  Preparation  and  coordination  of  the 
material  contained  herein  have  extended  over  the  past  six  years  with  contributions  from  the  various 
organizations  in  IRIG.  The  material  is  not  meant,  or  expected,  to  yield  complete  agreement  as  to  the 
relative  merits  or  importance  of  these  procedures.  However,  this  document  is  a  step  in  the  direction 
toward  the  eventual  establishment  of  IRIG  “guide  lines”  for  recommended  techniques  for  determining 
and  presenting  the  accuracy  estimates  of  data  collected  from  instrumentation  systems  in  support  of 
operational  testing.  In  addition,  this  volume  will  provide  a  source  of  documentation  on  current  and 
available  procedures  as  well  as  definitions  concerning  error,  accuracy,  and  precision. 
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PREFACE 


The  Inter-Range  Instrumentation  Group  has  had  as  one  of  its  goals  the  establishment  of  an  IRIG 
standard  for  techniques  in  determining  the  accuracy  estimates  for  instrumentation  systems  and 
methods  of  presenting  these  estimates.  The  IRIG  Steering  Committee  assigned  the  task  of 
accomplishing  this  goal  to  the  Data  Reduction  and  Computing  Working  Group  (DR&CWG). 

During  the  past  several  years  DR&CWG  has  endeavored  to  investigate  the  many  aspects  of  this 
task.  Work  has  proceeded  slowly  because  the  formulation  of  standard  and  acceptable  methods  for 
determining  accuracies  is  not  easy.  Furthermore,  it  is  difficult  to  obtain  agreement  on  the  means  of 
presenting  the  accuracy  estimates.  This  work  has  required  extensive  communication  and  coordination 
on  techniques  in  use  at  the  various  member  ranges.  As  a  result  of  this  continued  exchange  of 
information,  a  common  language  and  increased  understanding  have  evolved.  This  evolution  has  been 
slow,  but  it  now  appears  that  sufficient  agreement  between  member  ranges  and  adequate  technical 
maturity  have  developed. 

In  an  effort  to  get  some  results  completed  in  print,  the  DR&CWG  decided  to  break  the  task  into 
two  parts.  The  fust  portion,  dealing  with  the  positional  data  accuracies,  appeared  as  IRIG  Document 
Number  104-62.  The  second  portion,  scheduled  to  deal  with  velocity  and  acceleration  data  accuracies, 
was  written  in  1963  and  coordination  was  undertaken.  However,  at  the  21st  meeting  of  DR&CWG  it 
was  concluded  that  there  was  much  overlap  in  the  two  documents  and  that  it  would  be  desirable  to 
consolidate  them  since  104-62  was  then  at  least  three  years  old.  The  resulting  document,  103-64,  was 
greatly  amplified  and  contained  many  new  topics. 

This  new  document  basically  contains  the  material  in  103-64,  with  amplified  material,  and 
includes  some  of  the  newer  techniques  developed  and  used  during  the  past  few  years. 
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2.0 


INTRODUCTION 


To  assess  the  accuracy  of  the  final  result  of  processing  data  from  range  instrumentation  is 
extremely  difficult.  To  obtain  agreement  among  the  various  member  ranges  on  a  standard  for 
techniques  is  even  more  difficult.  The  major  reasons  for  these  difficulties  stem  from  the  variety  and 
diversity  of: 

1)  Various  types  and  amounts  of  instrumentation  used 

2)  Application  of  these  instruments 

3)  Operational  techniques  employed 

4)  Data  handling  procedures 

5)  Mathematical  and  computational  techniques 

6)  Types  of  requirements  for  reduced  data 

Although  the  member  ranges  have  a  wide  variation  in  instrumentation  and  operational 
techniques,  it  is  generally  agreed  that  the  measurement,  X(,  collected  at  an  instrument  at  time,  t,  is 
composed  of  a  “true  value”  or  “signal,”  and  an  error  component,  e{.  Thus 

X,  =  + 


or 


(2.0.1) 


et  Xt ' ^t 

The  common  problem  in  data  reduction  (and  error  estimation)  is  to  devise  mathematical  models  and 
corresponding  numerical  analysis  techniques  which  will  effectively  separate  the  signal  from  the  error  in 
some  “optimum”  manner. 

In  discussions  and  conferences  pertaining  to  error  estimation,  much  has  been  said  about  the  “true 
value”  of  a  parameter.  First,  it  is  emphasized  that  the  true  value  of  a  parameter  being  measured  is 
unknown.  If  the  true  value  were  known  there  would  be  no  need  for  further  discussion  of  errors  or  for 
taking  any  measurements  in  the  first  place.  As  Cassius  J.  Keyser  has  remarked,  “Absolute  certainty  is  a 
privilege  of  uneducated  minds  --  and  fanatics.”  Even  though  the  true  value  cannot  be  known  exactly,  it 
is  still  a  very  useful  concept  which  is  used  in  the  construction  of  mathematical  models  to  represent  and 
estimate  the  error  in  a  measurement  or  set  of  measurements. 

Since  the  true  value  of  the  signal  is  unknown,  the  approach  is  to  devise  mathematical  models  and 
experimental  design  techniques  to  give  “good”  or  “best”  estimates  of  the  signal.  In  data  reduction 
terminology,  these  estimates  are  often  referred  to-  as  “computed  values”  or  as  “standards  for 
comparison”  in  error  analysis.  The  computed  value,  /i£,  is  a  numerical  function  of  the  measurements 
X(  or  may  be  a  simultaneous  measurement  from  a  more  accurate  instrumentation  system  or  several 
systems.  The  estimate  of  the  error  is 
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et*et  =  V'V 


(2.0.2) 


The  quantity  e£  is  called  a  residual,  and  applied  error  analysis  utilizes  this  quantity  as  its  basic 
input  in  estimating  error  characteristics. 

It  is  important  to  note  that  the  measurement  Xt  is  a  random  variable  and  since  and  et  are 
numerical  functions  of  X£,  they  are  also  random  variables.  Each  of  these  quantities  are  functions  of 
time  and  thus  make  up  a  realization  or  sample  from  a  stochastic  or  random  process.  If  the  observations 
are  taken  at  discrete  times,  they  form  a  sequence  of  random  variables  or  a  stochastic  sequence  and  are 
commonly  referred  to  as  a  discrete  time  series.  This  concept  readily  suggests  that  probability  and 
statistics  play  an  important  role  in  error  estimation,  as  well  as  in  obtaining  parameter  estimates  from 
the  data  reduction  process. 

The  concept  of  accuracy  with  respect  to  measured  and/or  reduced  data  is  closely  related  to  the 
error  therein,  but  is  not  identical  to  it.  Certainly  “small  error”  implies  high  accuracy  or  “accurate” 
data,  and  “large  error”  corresponds  to  low  accuracy  or  inaccurate  data.  In  general  it  is  sufficient  to 
state  that  accuracy  is  some  function  of  the  error  distribution,  and  in  the  final  analysis,  accuracy  itself 
must  be  estimated  and  based  on  quantities  which  are  approximated  numerically,  assumed,  or  based  on 
statistical  estimation  of  parameters  which  characterize  the  error  distribution. 
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3.0 


ERROR  CLASSIFICATION,  ACCURACY,  AND  PRECISION 


In  accordance  with  Section  2.0  the  error  is  estimated  by  the  residual  ef  and  the  accuracy  of  the 
measured  and/or  reduced  data  is  estimated  by  some  function  of  the  residual.  This  function  usually 
involves  the  expected  value  of  the  error.  The  expected  errors  are  generally  classified  into  two  basic 
types,  the  random  (noise)  errors  and  the  systematic  errors.  The  theoretical  error  model  is 

et*St+Nt  (3.0.1) 

where 

S£  =  systematic  error,  and 

N£  =  random  or  noise  error. 

The  random  error  is  nondeterministic.  This  means  it  is  not  described  by  an  analytic  function  but 
must  be  characterized  in  terms  of  its  probability  distribution  function.  The  systematic  error  is 
generally  considered  to  be  deterministic  and  can  be  represented  with  an  analytic  model.  If  the  error 
distribution  is  thought  of  in  terms  of  variation,  then  the  systematic  error  comprises  the  “explained 
variation”  while  the  random  error  corresponds  to  the  “unexplained  variation.” 

Another  type  of  error  which  is  not  explicitly  covered  by  the  model  (3.0.1),  but  which  cannot  be 
overlooked  is  “gross  error.”  This  type  of  error  is  usually  due  to  some  type  of  malfunction  in  the 
measurement  process  and  is  not  a  usable  or  valid  measurement.  For  this  reason  data  with  this  type  of 
error  is  referred  to  as  “bad  data,”  “wild  points,”  or  more  formally  as  “outlyer  data.”  This  type  of 
error  does  not  belong  to  the  statistical  population  described  by  the  probability  distribution  function 
of  Nt  and  hence  must  be  removed  or  edited  out  so  it  will  not  bias  or  corrupt  the  data  reduction 
process.  Many  different  types  of  editing  techniques  are  used  on  the  various  ranges  and  are  not 
discussed  in  this  report.  However,  it  is  emphasized  that  if  outlyers  are  ignored  in  the  data  reduction 
process  their  effect  will  corrupt  reduced  data  from  neighboring  valid  measurements  and  render  the 
otherwise  valid  data  unusable. 

It  is  important  to  note  that  it  may  often  be  difficult  to  distinguish  random  from  systematic  error 
in  a  given  set  of  measurements.  This  may  be  due  to  the  fact  that  certain  constraint  conditions  exist 
(geometry  limitations,  availability  and  capability  of  instrumentation,  cost,  safety,  etc.)  and  no  feasible 
test  can  be  designed  which  will  allow  a  reduction  process  to  give  adequate  residuals;  or  it  may  be  that 
St  is  very  complicated  and  an  adequate  error  model  is  not  available;  or  it  may  be  that  the  variance  of 
the  random  error  is  large  and  the  ratio  of  systematic  error  to  noise  error  is  small  (error  signal  to  noise 
ratio),  which  corrupts  the  estimates  made  of  St.  Nevertheless,  in  concept,  and  often  in  practice,  it  is 
possible  to  treat  random  errors  statistically  while  the  systematic  errors  may  be  estimated  in  a 
deterministic  manner  and  corresponding  corrections  made  to  the  observed  data. 

The  combination  of  systematic  and  random  error  is  referred  to  as  total  error,  e£,  and  in  this 
document  the  accuracy  of  a  measurement  is  defined  as  the  root  mean  square  of  the  distribution  of  the 
total  error.  In  statistical  terms 


(3.0.2) 

Accuracy  «.  i/e  [(total  error)4]  ■ 
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where  E  (•)  is  the  expected  value  of  the  quantity  in  brackets. 


If  the  systematic  error  and  random  error  are  independent,  then 
Accuracy  *  A  [(systematic  error)2]  +  E  [(random  error)2]. 


(3.0.3) 


Another  important  term  which  should  be  distinguished  from  accuracy  is  “precision.”  There  are 
probably  few  words  as  loosely  used  by  scientific  personnel  as  precision  and  accuracy.  Accuracy  can  be 
described  as  “closeness  to  the  truth”  while  precision  is  the  “closeness  together”  or  internal  consistency 
of  a  set  of  measurements  (see  Figure  1 ).  From  these  descriptions  we  define  precision  as 


Precision  ~  /t  [(random  error)2]. 


(3.0.4) 


It  is  seen  that  accuracy  requires  precision  but  precision  does  not  necessarily  imply  accuracy.  If 
the  systematic  error  >  random  error  then  in  this  special  case 


Accuracy  ~  / fc  ((systematic  error)2]. 


(3.0.5) 


It  is  recognized  that  among  the  ranges  there  is  not  necessarily  unanimous  agreement  regarding  the 
classification  of  errors  and  the  corresponding  concepts  of  precision  and  accuracy.  However,  in  the 
interests  of  practicability  the  errors  are  classified  as  random  and  systematic  with  the  important 
distinction  that  the  term  “bias  error”  is  considered  to  be  systematic  error  and  “constant  bias  error”  is 
considered  to  be  the  constant  component  of  systematic  error. 

At  times  it  is  convenient  to  characterize  and  analyze  error  in  the  frequency  domain.  When  this  is 
done  for  trajectory  type  data  the  error  classification  may  be  grouped  into  two  categories.  The  high 
frequency  category  may  correspond  in  a  sense  to  random  noise  errors  while  the  low  frequency  end  of 
the  spectrum  may  correspond  to  systematic  errors  (although  not  necessarily).  The  systematic  error  has 
a  special  case  corresponding  to  zero  frequency  which  is  constant  bias  error.  The  frequency  spectrum  is 
a  useful  tool  in  the  analysis  of  error  and  is  discussed  in  more  detail  in  Section  1 1 .0. 

To  effectively  estimate  systematic  errors  with  the  idea  of  predicting  what  accuracy  can  be 
expected  on  future  tests,  reliance  must  be  made  on  experience  from  past  error  analyses  concerning  the 
variation  of  systematic  errors  within  a  test  and  from  test  to  test.  Viewed  in  this  manner,  the  systematic 
errors  are  random  variables  which  may  be  characterized  by  a  probability  distribution  function.  An 
example  of  how  the  constant  bias  estimate  from  a  residual  distribution  of  an  instrumentation  system 
such  as  a  radar  will  vary  from  test  to  test  is  illustrated  in  Figure  J.a. 

The  variation  of  errors  within  and  between  tests,  and  between  instrumentation  sites  leads  one  to 
look  for  techniques  which  will  determine  whether  or  not  significant  differences  exist  in  the  error  from 
variable  factors  and  to  estimate  what  effects  changes  in  these  factors  have  on  the  error.  Techniques 
relating  to  this  type  of  error  analysis  deal  with  the  analysis  of  variance  and  are  discussed  in  Section 
6.3. 
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It  is  recognized  that  there  are  other  descriptions  which  may  be  applied  to  errors.  Correlation  of 
errors  in  all  of  its  detail  may  not  be  fully  understood  at  this  time.  The  various  interdependencies  of 
data  and  the  added  complications  which  arise  from  auto,  serial,  and/or  time  series  correlation  makes  it 
unwise  to  attempt  a  full  classification  at  this  time;  however,  the  importance  of  correlation  between 
and  within  errors  should  not  be  ignored  and  a  general  discussion  of  correlation  and  its  effects  on  error 
estimation  is  given  in  Sections  3.3  and  4.7. 
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T 


1  -  Low  Accuracy, 

Low  Precision 


2  -  Low  Accuracy, 

High  Precision 


3  -  Low  Precision, 


No  Systematic  Error 


4  -  High  Accuracy, 

High  Precision 


5  -  Systematic  Error 

FIGURE  1.  The  concepts  of  accuracy,  precision  and  error  are  illustrated. 
Consider  T  to  be  a  target  with  the  dots  representing  the  results  from  several  * 
firings  at  the  target 
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tlmt  (secs )  -* 

FIGURE  l.a  The  systematic  error*  Z  in  azimuth  A,  as  measured  by  a  radar, 
have  variations  as  pictured  above  both  from  test  to  test  and  as  a  function 
of  time  for  a  particular  test. 


3.1  ERRORS  IN  BASIC  QUANTITIES 


Thus  far  only  errors  in  measurements  or  observations  have  been  discussed.  However,  a  great 
variety  of  error  sources  enter  into  the  accuracy  considerations  for  data  systems  Lied  in  tracking 
missiles.  At  this  time  and  within  the  present  state-of-the-art,  it  is  not  practical  or  even  possible  to 
accurately  estimate  the  error  in  position,  velocity,  and  acceleration  due  to  each  of  these  error  sources. 
In  an  attempt  to  meet  demands  for  more  and  more  stringent  accuracy  requirements,  great  effort  has 
been  directed  toward  the  improvement  of  instrumentation.  However,  it  must  be  realized  that  even 
with  perfect  instrumentation,  there  are  basic  theoretical  limitations  on  the  accuracy  which  may  be 
achieved.  These  limitations  may  introduce  errors  which  make  impossible  the  accuracy  which  is  being 
demanded  of  the  tracking  systems. 

Typical  errors  associated  with  some  basic  quantities  are  listed  below.  These  exist  quite  apart  from 
the  instrumentation  and  place  limitation  or  lower  bound  on  the  accuracy  which  may  be  attained.  It  is 
realized  that  there  may  be  no  unanimity  regarding  these  magnitudes,  or  for  that  matter,  their  types, 
i.e.,  probable,  absolute,  etc.,  but  they  will  be  given  as  points  of  departure. 


a) 

First  order,  Class  I  surveys  (distance) 

10  ppm  (parts/1000000) 

b) 

Velocity  of  light 

1-2  ppm 

c) 

Index  of  refraction 

25  ppm 

d) 

Changes  in  refraction  due  to  rapid 

fluctuations  in  atmosphere 

30-40  ppm 

e) 

Errors  in  direction  cosines  caused  by 

the  above  fluctuations 

10  ppm 

0 

Overall  ballistic  camera  accuracy 

5-15  ppm 

g) 

Errors  in  reading  ballistic  camera  film 
due  to  1°C  fluctuation  in  reading 

room  temperature 

5  ppm 

h) 

Star  catalog  error 

1  ppm 

i) 

Error  in  100  mile  base  line 

1-5  feet 

j) 

Error  in  origin  of  national  survey 

50-100  feet 

k) 

Differences  between  International 

Spheroid  and  Clarke  1866  Spheroid 

150-300  feet 

1) 

Undulations  of  the  geoid  (Aa) 

50-100  feet 

m) 

Capability  to  determine  relative  geoid 

heights  at  widely  separated  stations 

30  feet 

n) 

Voltage  level  (secondary  standard) 

100  ppm 

°) 

Mass 

.0001-100  ppm 

P) 

Timing  (atomic  standard) 

10” 12  seconds 

q) 

Timing  (laboratory  electronics) 

10"8  seconds 

r) 

Timing  (instrument  crystal) 

10"6  seconds 

s) 

Timing  correlation  between  separate  stations 

10*5  seconds 

Some  of  the  above  values  are  contained  in  Reference  (1). 
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It  is  interesting  to  observe  the  effects  of  the  quoted  survey  and  refraction  errors  when  propagated 
into  spatial  position  and  velocity  at  a  range  of  250  miles.  The  resulting  errors  due  to  these  sources 
alone  are: 

Position:  15  feet  in  horizontal  coordinates 
40  feet  in  vertical  coordinates 

Velocity:  0.1 5  foot/second  in  horizontal  coordinates 
0.40  foot/second  in  vertical  coordinates 

Obviously,  these  are  facts  which  must  be  considered  when  imposing  tracking  accuracy 
requirements  and  in  designing  instrumentation  systems  to  meet  them.  One  is  immediately  faced  with 
the  operation  of  the  law  of  diminishing  returns. 


3.2  CAUSES  OF  ERRORS  IN  INSTRUMENTATION  TRACKING  SYSTEMS 


Once  a  set  of  residuals  has  been  calculated  with  corresponding  estimates  of  the  random  and 
systematic  errors,  the  process  of  determining  the  cause  of  the  errors  is  often  extremely  difficult.  The 
diffculty  encountered  in  determining  the  cause  of  error  depends  on  several  factors  such  as  the  ability 
to  effectively  design  and  control  an  experiment,  knowledge  of  the  variables  which  effect  the  error, 
correlations  and  effects  between  variables,  correct  assumptions  concerning  error  model  and  probability 
distribution  functions. 

Some  causes  of  random  errors  in  measurements  may  be  tropospheric  and  ionospheric 
scintillations,  electronic  noise,  instrumentation  wear,  mechanical  play,  granularity,  or  tolerances  in  the 
measurement  capabilities  of  the  instrument,  errors  in  computing  and  timing,  multipath,  scintillation 
due  to  echo  or  skin  track,  Film  reading,  and  so  on. 

Systematic  errors  are  caused  by  physical  limitations  in  achieving  a  true  geodetic  survey, 
electronic,  optical  and  mechanical  misalignment,  servo  lag,  phase  and  frequency  drift,  frequency  and 
timing  bias,  encoder  nonlinearity,  antenna  droop,  mislevel,  lens  distortion,  beacon  delay, 
miscollimation,  dial  eccentricity,  and  nonorthogonality  of  azimuth  and  elevation  axes.  One  of  the 
most  important  types  of  systematic  error  is  usually  considered  to  be  constant  bias  or  zero  set  error. 
This  error  is  very  important  because  it  is  often  the  most  significant  component  of  the  systematic  error 
St.  In  addition,  it  is  the  easiest  component  to  estimate  and  to  correct  for;  in  fact,  under  actual  flight 
conditions,  the  dynamics  and  geometry  are  often  of  such  nature  that  the  constant  bias  error  is  the 
only  component  of  the  systematic  error  which  can  be  effectively  estimated.  Because  of  its  constant 
nature,  attempts  to  correct  the  instrumentation  for  this  type  of  error  can  be  made  using  static 
calibration  tests  before  and  after  the  actual  flight  tracking  operation.  Although  constant  bias  error  in 
the  raw  measurements  can  be  considered  practically  constant  for  the  duration  of  a  single  test,  its 
propagated  effects  on  position  error  may  vary  considerably  during  a  single  test  because  of  sensitivity 
to  tracking  geometry.  Neighboring  spatial  points  in  the  time  and  space  domains  will  in  general  have 
nearly  constant  error  due  to  a  Fixed  bias  in  misalignment.  However,  propagated  position  error  would 
be  significantly  different  at  t+10  seconds  (close  range)  from  what  it  would  be  at  t+150  seconds  (slant 
range  of  several  hundred  miles).  However,  position  errors  due  to  the  misalignment  in  the  period  t+140 
to  t+160  seconds  would  be  expected  to  remain  almost  constant  depending  on  geometry  changes  in 
that  particular  interval. 


3.3  CORRELATION  OF  RANDOM  ERROR 


The  random  errors  Nj,...,Nn  corresponding  to  measurements  Xi,...,Xn  may  be  interrelated  so 
that  the  error  in  the  jth  measurement  may  be  dependent  upon  the  errors  in  X:  ,  ,Xj_2,...Xj_p,  p<j.  If 
this  relationship  exists,  the  errors  are  dependent  and  said  t»>  be  correlated.  In  tne  case  of  a  time  series 
this  interdependence  is  called  the  autocorrelation  and  is  estimated  by  means  of  computing  the  serial 
correlation. 


From  a  mathematical  standpoint  two  random  variables  X  and  Y  are  independent  if  their  joint 
probability  density 

F(x,y)  =  f(x)g(y) 

where  f(x)  and  g(y)  are  the  probability  density  functions  corresponding  to  X  and  Y  respectively.  Two 
random  variables  are  uncorrelated  if  E(XY)=E(X)E(Y).  It  can  be  shown  that  independence  implies  two 
random  variables  are  uncorrelated,  but  not  vice  versa.  The  correlation  coefficient  between  the  random 
variable  X  and  Y  is  defined  as 


E  [  (X  -  px)(Y  -  uY)] 


(3.3.0) 


[E(X  -  yx)z  E(Y  -  yY)2P  uXuY 

From  relation  (3.3.0)  we  see  the  correlation  coefficient  is  the  ratio  of  the  covariance  to  the  product  of 

the  standard  deviations  of  X  and  Y.  If  a  sample  X,  ,...,Xn,  Y, . Yn  is  taken,  the  correlation 

coefficient,  pxy  is  estimated  by  „  V  «.  -  XHY.  -  ?) 


cuciiicicni,  pxy,  t.  cMimaicu  oy  (N  .  j)  J  (Xj  -  X)  (Yj  -  Y) 

- - - - ; - i-  (3.3.1) 

xy  N  N  'j'i 

N  l  -  X)2  l  (Y1  -  Y)2 

'“i-1  i-1 

lfX|,X2,..  Xn  is  a  realization  the  autocorrelation  p(k)  at  lag  k  is  estimated  by  the  serial 
correlation 


? TTk  l  «i-*HX1+k-x> 

_ >1 _ 

r  N  N-k 

E  <xi  "  x^2  N  -  k  -  1  ^  (Xi 

i-1  i-1 


(3.3.2) 


k  -  0,1.2 . M<N 
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(3.3.2)  is  analogous  to  (3.3.0)  in  that  the  autocorrelation  function  is  the  ratio  of  the  autocovariance  to 
the  standard  deviations  of  Xj  and  Xj+j^.  If  we  let  the  autocovariance  estimate  at  lag  k  be  R(k)  and  if 
the  stochastic  process  is  stationary,  then  the  correlation  coefficient  can  be  replaced  by  the  relation 

n=Ml  (3.3.3) 

R(0) 

(3.3.3)  is  called  the  normalized  auto-covariance  or  the  circular  correlation  coefficient.  The  graph  of  r^ 
is  sometimes  called  a  correlogram  (see  Figure  3). 

The  autocorrelation  between  errors  is  especially  known  to  be  present  if  the  measurements  are 
from  an  instrument  whose  measurements  are  taken  with  a  small  time  increment  between  them.  Most 
types  of  range  instrumentation  are  of  this  type  and  as  a  result  the  autocorrelation  of  the  errors  occurs 
frequently. 

The  effects  of  high  correlation  between  errors  can  lead  one  to  trouble  in  both  data  reduction  and 
error  analysis.  If  the  correlation  among  errors  is  high  and  does  not  “die  out”  or  become  small  in  a 
short  time  the  observed  data  will  deceptively  appear  to  be  smoother  (see  Figure  2),  or  to  contain  what 
appears  to  be  low  frequency  oscillations  and  trends  which  could  very  easily  be  falsely  assumed  to  be 
“signal”  or  valid  information,  and  could  lead  one  to  false  conclusions  about  the  proper  technique  for 
editing  and  smoothing  of  the  data.  In  addition,  it  can  be  shown  that  if  high  correlation  exists,  then 
estimates  of  the  variance  of  the  error  distribution  will  be  biased  (see  Section  4.7). 

Several  methods  are  used  for  modeling  autocorrelated  data.  The  general  idea  is  to  express  the 
random  error  as  a  stochastic  function  which  depends  on  parameters  which  are  random  variables.  The 
most  common  is  a  pth  order  stochastic  difference  equation  or  pth  order  Markov  process.  The  first 
order  Markov  process  gives  a  probabilistic  model  of  the  form 

Ni  =  pNi.,+7ji  (3.3.4) 

where  j?j  is  an  independent  random  variable  0  <  p  <  1.  The  serial  correlation  coefficient  corresponding 
to  this  process  can  be  shown  to  be  exponential,  i.e.,  the  correlation  between  the  errors  in  the  ith  and 
(i+k)th  measurement  is  p^,  where  p  is  the  correlation  coefficient  between  adjacent  measurements. 
Figure  2  shows  independent  noise  (p= 0)  and  correlated  noise  corresponding  to  a  first  order  Markov 
process  with  p=. 3,  p=. 6,  and  p=. 99,  respectively.  Inspection  of  the  graphs  of  the  errors  in  Figure  2 
shows  that  as  the  correlation  becomes  high,  the  data  appears  smoother  and  systematic  trends  appear  fn 
the  error  profile  which  can  be  very  misleading. 

The  second  order  Markov  process  is 

Ni  =  a»Ni-,  +a2Ni-2  +,?i 

where  |atf  <1,  |a2|  <1,  and  rjj  is  an  independent  random  variable.  In  this  model  the 
autocorrelation  damps  out  exponentially  but  also  contains  a  sinusoidal  oscillation.  An  example  of 
FPS-16  error  estimates  which  follow  this  type  of  autocorrelation  is  given  in  Figure  3.  Note  that  on  this 
figure  the  autocorrelation  is  low  at  a  lag  of  0.5  seconds  between  points  (i.e.,  the  correlation  between 
points  separated  by  0.5  seconds  is  low). 
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FIGURE  2  Example  of  uncorrelated  noise  (RHO  =  0)  and  low  correlated 
noise  (RHO  =  .3)  to  very  high  correlated  noise  (RHO  *  .99). 
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FIGURE  3  CORRELOGRAM  OF  AN/FP5 -  16  RANGE,  AZIMUTH 
ANO  ELEVATION  RESIDUAL  DATA  USING  LAG 
OF  .5  SECONDS. 


4.0 


METHODS  USED  TO  ESTIMATE  RANDOM  ERRORS 


From  the  previous  discussion  it  is  seen  that  random  errors  are  estimated  in  terms  of  parameters 
which  describe  or  characterize  their  behavior  in  terms  of  a  probability  distribution  function.  The 
parameters  estimated  are  almost  always  in  terms  of  the  first  and  second  sample  moments 

corresponding  to  the  residual  distribution.  If  e, ,  e2 . en  Is  a  set  °f  errors,  then  the  variance  of  the 

random  error  is  expressed  as 


o*  -  E(eJ)  -  [E(£i)]2,  (4.0.1) 


and  if  e1 ,  e2 . en  is  a  set  of  n  residuals  from  a  stationary  time  series,  which  estimates  this  set,  then  the 

random  error  is  estimated  by 


6 


e 


(4.0.2) 


where 


n 


i-1 


Relation  (4.0.2)  is  the  sample  standard  deviation  of  the  residual  distribution.  In  general,  the  trajectory 
measurements  themselves  are  changing  and  do  not  form  a  stationary  process.  However,  the  random 
error  of  the  measurements  is  usually  assumed  to  be  stationary  random  process.  This  means  that 
effective  methods  for  removing  nonstationary  effects  must  be  devised.  Thus,  the  problem  in  estimating 
the  random  error  is  to  be  able  to  separate  the  nonstationary  signal  from  the  error,  obtain  a  set  of 
residuals,  then  separate  the  systematic  error  from  the  random  error.  This  process  is  dependent  on 
many  factors  including  such  elements  as  computational  techniques,  autocorrelation  and  cross 
correlation  and  assumptions  concerning  the  analytic  form  of  the  systematic  error  model. 

It  is  especially  desirable  to  have  a  good  estimate  of  the  random  errors  for  a  particular  system  on  a 
particular  test  as  a  function  of  time.  This  information  is  needed  to  settle  such  important  questions  as 
whether  or  not  the  successive  errors  in  the  measured  quantity  are  correlated,  what  the  optimum 
smoothing  functions  are,  and  what  the  best  methods  are  for  estimating  velocity  and  acceleration  data 
from  position  data. 

Assuming  that  the  autocorrelation  and  the  cross  correlation  are  zero,  several  methods  of 
estimating  random  errors  may  be  considered,  and  are  given  in  Sections  4.1  through  4.6. 


4.1  GRAPHICAL  TECHNIQUES 


In  some  instances  graphical  profiles  of  the  residual  distribution  versus  time  are  available  for 
analysis.  This  is  especially  true  of  electronic  measuring  equipment  where  the  servo  feedback  voltages 
are  recorded  graphically  in  analog  form.  In  cases  such  as  this,  the  residual  profile  can  be  inspected  to 
see  if  there  are  systematic  trends,  etc.  If  the  residual  data  is  free  of  trends  due  to  systematic  error  and 
correlation,  one  can  select  a  fixed  or  variable  sample  rate  and  use  relation  (4.0.2)  in  estimating  the 
random  error. 

If  the  residual  is  from  a  population  characterized  by  the  normal  distribution  then  a  much  simpler 
method  is  available  and  is  based  on  peak-to-peak  or  the  sample  range  over  a  selected  group  or  block  of 
residuals.  In  selecting  groups,  we  take  k  groups  each  containing  n  points.  The  ranges  Rj  which  are 
peak-to-peak  values  for  each  jth  group  are  estimated  by 

Rj  =  max  ey  -  min  ey  (4.1.1) 

The  average  range  R  of  the  sample  ranges  !s  rl  »n  computed,  and  the  standard  deviation  is 
estimated  from  the  relation 

a  =  anR  (4.1.2) 

where  an  is  the  ratio  o/E(R)  for  the  standard  normal  distribution  and  is  available  in  tables  of  the 
distribution  of  the  standardized  range  for  a  normal  population  (see  Reference  23  for  tables). 

Since  the  peak-to-peak  technique  is  relative  to  the  “envelope”  of  the  residual  distribution  for 
each  group  or  subinterval,  the  size  of  the  subinterval  can  be  selected  so  that  trends  can  be  removed  if 
they  exist  in  the  data. 

The  graphical  techniques  can  be  very  useful  for  “quick  look”  analysis  but  have  a  drawback  in 
that  they  usually  involve  much  manual  analysis.  When  using  this  technique  one  must  be  careful  that 
“wild  points”  are  omitted  since  the  range  technique  is  not  as  stable  as  estimating  a  by  relation  (4.0.2). 
In  addition,  one  must  be  certain  that  the  serial  correlation  between  sampled  points  is  not  high. 
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4.2 


LEAST  SQUARES  POINT  ESTIMATES  OF  RANDOM  ERRORS 


This  method  can  be  most  easily  used  when  a  point-by-point  least  squares  fit  to  a  redundant  data 
set  is  being  used  for  reduction  of  the  data. 

Examination  of  the  discrepancies  between  the  solutions  and  the  observed  measurements  will 
yield  information  regarding  the  distribution  of  errors  in  the  measurements.  By  combining  many  of 
these  residuals  one  may  estimate  their  standard  error.  A  comparison  of  this  estimate  with  that 
predicted  by  the  a  priori  statistics  of  the  adjustment  gives  an  approximation  to  the  error  in  this  a  priori 
variance  data.  By  repeating  this  procedure  using  new  estimates  for  the  a  priori  variance  data,  an 
improved  figure  for  the  data  error  may  be  found. 

The  sample  standard  deviation  of  the  space  position  residual  at  each  ith  time  point  is 


9 


i 


(4.2.1) 


where 

ey  is  the  residual  distance  at  the  ith  time  for  the  jth  instrument 
k  =  number  of  measured  parameters  (R,  A,  E,  etc.) 
n  =  number  of  instruments 

3  =  number  of  parameters  being  estimated  {X,  Y,  Z). 

The  term  kn-3  represents  the  number  of  degrees  of  freedom  available  in  the  estimate  of  the  variance. 
In  the  case  of  an  optical-only  solution  the  degrees  of  freedom  would  be  2n-3,  and  3n-3  for  a  solution 
using  n  radars.  In  the  case  where  nj  radars  and  n2  optical  instruments  were  used,  we  would  have 
3nI  +2n2-3  degrees  of  freedom. 


Position  error  estimates  in  the  required  coordinate  system  are  almost  always  correlated  and  as  a 
result,  the  position  error  distribution  is  characterized  by  the  estimate  of  the  variance-covariance  or 
“dispersion”  matrix.  The  variance-covariance  matrix  estimate  for  the  ith  time  point  is 
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where  N'1  is  the  inverse  coefficient  matrix  corresponding  to  the  normal  equation  obtained  by 
computing  the  partial  derivatives  in  the  least  square  process. 

Care  must  be  taken  in  applying  this  method,  because  it  is  capable  of  detecting  systematic  or  bias 
error  as  well  as  random  error.  An  examination  of  the  time-dependent  behavior  of  the  residuals  should 
be  made  to  determine  whether  or  not  apparent  serial  correlations  or  trends  among  these  errors  are  due 
to  systematic  causes. 

If  corresponding  residuals  from  point  to  point  display  “trends,”  residual  systematic  errors  are 
indicated  which  often  may  be  estimated  and  removed  by  techniques  such  as  a  “best  estimate  of 
trajectory”  (BET)  solution,  a  multiple  regression  analysis  or  numerically  filtering  over  a  series  of 
points.  After  the  trends  or  systematic  error  estimates  have  been  removed,  estimates  of  random  error 
may  be  made  utilizing  the  new  residuals  which  correspond  to  the  “unexplained”  or  random  variations 
about  the  trend. 

In  cases  where  all  observations  have  equal  weight,  the  standard  error  in  the  observations  may  be 
estimated  by  (4.2.1).  For  the  weighted  case  the  quadratic  form  of  the  dispersion  matrix  of  the 
residuals  mav  be  used  to  obtain  an  estimate  of  the  unit  variance.  That  is  a^=e.’W',e  where  e  is  a  column 
vector,  W  a  weight  matrix.  In  this  case  the  dispersion  matrix  would  be 

Dj  =  (T-  Wj*  Tj)*1  of  (4.2.3) 


In  (4.2.3)  T  is  the  “design”  matrix,  T’  is  its  transpose  and  the  coefficient  matrix  of  the  normal 
equations  is  N=TT.  The  matrix  Wj  is  the  weight  matrix.  The  consistency  of  this  estimate  may  be 
established  by  comparing  its  diagonal  terms  with  a  value  of  the  variances  known  a  priori. 


4.3  MULTI-INSTRUMENT  ESTIMATE  OF  RANDOM  ERROR 


This  method,  sometimes  called  Simon-Grubbs,  estimates  random  error  using  over  determination 
from  additional  independent  systems.  In  these  computations  (simplified)  the  mathematical  model  for 
three  measurements  is 


yji  =xji  +eji, 

yj2  =xj2  +ej2, 
yi3  =x}3  +ej3_ 


where  y-  is  the  measured  value  of  the  ith  characteristic  xy  with  associated  measurement  error  eij for 
the  jth  instrumentation  system.  Form  the  three  differences: 


Ai*2  =  ej,  -ej2 
^1.3  =  eJi  "  ej3 
A2,3  =  eJ2  •  ei3* 

The  bias  error  is  tacitly  assumed  to  be  some  fixed  value;  the  variance  estimates  for  the  above 
difference  equations  are  computed  about  the  mean  difference.  If  the  bias  is  changing,  it  adds 
components  of  variance  estimates  given  below: 

s?,2  =  s2,  +  s22 

s?,3  = 

Sj,3  =  Sgl  +  Sg3 


The  above  equations  are  solved  for  s2  j ,  s2  2 ,  and  s2  3 ,  the  respective  error  variance  estimates  for 
instrumentation  systems  1,  2,  and  3.  The  variance  of  the  estimated  error  variance  for  instrumentation 
System  One  is  given  by: 


2s 


Var 


ei 


(n  -  1)  T  (n  -  1) 


(s2  s2 
ei  e2 


s2  a2  + 
ei  es 


i2 

e2 


es 


For  measurement  Systems  Two  and  Three,  similar  expressions  with  the  subscripts  permuted  are 
appropriate. 


4.4  VARIATE  DIFFERENCE  METHOD 


The  variate  difference  method  is  a  technique  which  uses  the  pth  order  successive  difference  to 
“detrend”  data  and  estimate  the  variance  of  the  random  error  about  the  trend.  It  can  be  shown  that 
the  coefficients  aj  of  the  linear  combinations  for  the  pth  difference  follow  the  binomial  expansion  law 
and  thus  the  sum  of  the  a?  for  the  pth  difference  is 


P+1 


l 

i-1 


_ £2e1! _ 

p!(2p  -  p)! 


(4.4.1) 


Since  successive  differences  are  an  analog  to  the  time  derivative  it  is  obvious  that  constant  linear  and 
quadratic  trends  will  be  removed  from  the  measured  data  by  a  first,  second  and  third  order 
differencing  scheme.  In  actual  practice  the  third  difference  is  usually  sufficient  and  in  almost  all  cases 
the  third,  fourth  and  fifth  difference  estimates  give  equivalent  results. 

If  X,  ,...Xn  is  a  set  of  measurements  or  sample  from  a  set  of  independent  random  variables  with 
mean  jux  and  variance  a*  then  the  linear  combination 

n 

Y.  ■  l  VX1  (4.4.2) 

i-1 


is  a  random  variable  with  mean 


f  n  1 


(4.4.3) 


and  variance 


(4.4.4) 
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When  taking  successive  differences,  the  resulting  difference  turns  out  to  be  a  linear  combination 
of  the  original  random  variable.  Each  difference  gives  a  new  random  variable  which  can  be  expressed  as 
a  sum  of  the  original  measurements  {  Xj  }.  For  example,  the  first  difference  is 


Y  -  A1 
li  x. 


Xi  ’  Xi-,  • 


*1  ■  1.  *2 


1 


For  the  second  and  third  differences  we  get 

*21  *  ^  *  *1+1  '  2X1  +  Xi-,  >  *.  *  1.  »2  -  -  2,  -  1  , 

Y>1  “  4xt  ”  Xi+2  '  3X1+1  +  3Xi  '  Xl-1  .  *1  “  1.  »2  "  -3.  *3  -  3,  -  -1 

and  so  on. 


From  (4.4.4)  and  (4.4.1)  the  equation  for  the  variance  of  the  pth  difference  is 


(4.4.5) 


Now  the  estimate  of  o2 


is 


so  that 


where 


e2 

X 


(pp)  fa  -  p> 


A£,=  pth  successive  difference  in  Xj 
n  =  number  of  data  points  in  sample 
p  =  order  of  the  successive  difference 


(4.4.6) 


(4.4.7) 


The  variate  difference  method  is  widely  used  since  it  can  be  applied  to  data  when  no  knowledge 
concerning  the  trend  in  the  data  is  available.  However,  caution  must  be  used  because  of  the 
assumptions  which  were  made  concerning  independence  of  the  errors.  It  should  be  emphasized  that  if 
the  error  is  highly  correlated,  then  the  error  variance  estimates  will  be  biased,  and  in  fact  usually 
underestimated  (see  Section  4.7).  In  addition  to  the  correlated  effects  on  the  variate  difference 
method  the  technique  is  also  highly  unstable  if  “wild  points”  or  outlyers  are  sampled.  One  method  of 
compensating  for  autocorrelation  in  the  data  is  to  compute  the  differences  using  a  larger  time  spacing 
At  between  the  measured  Xj’s.  However,  one  must  be  careful  that  when  large  time  increments  are 
used,  the  effects  of  aliasing  do  not  bias  the  variance  estimates.  A  more  satisfactory  technique  is  to 
estimate  the  autocorrelation  function  and  use  this  to  correct  the  biasing  effects  on  (4.4.7).  The 
amount  of  biasing  caused  by  autocorrelation  is  discussed  in  Section  4.7. 

Since  the  variate  difference  method  is  based  on  differencing  it  is  a  useful  tool  in  examining  the 
frequency  content  of  pth  order  differences  obtained.  If  At  is  the  time  interval  between  successive  data 
points,  the  maximum  distinquishable  frequency  is  (l/2At)  cycles  per  second  and  is  called  the  Nyquist 
frequency.  Higher  order  differences  around  the  order  of  4,  5  or  6  effectively  examines  the  error 
contributions  of  frequencies  between  0.82/2At  to  the  Nyquist  frequency.  This  is  due  to  the  fact  that 
frequencies  below  0.82/2At  have  been  removed  by  differencing.  If  we  can  assume  the  noise  is  uniform 
over  all  frequencies,  we  may  estimate  the  component  of  the  variance  in  the  interval  .82/2At  to  l/2At 
by  dividing  the  variance  obtained  using  the  variate  difference  method  by  four.  The  frequency  range  of 
0.82/2At  to  l/2At  is  based  on  the  frequency  response  of  the  variate  difference  method.  For  a  fourth 
order  difference  the  weights  are  (1,  -4,  6,  -4,  1).  If  we  take  the  Fourier  transform  of  these  weights  we 
obtain  the  frequency  response  as  given  in  Figure  4. 


Figure  4  -  Frequency  response  for  fourth  order  variate  difference  method  with 
lag  At  between  successive  points. 
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4.5  THE  LEAST  SQUARE  CURVE  FITTING  METHOD 


In  the  least  square  method  a  polynomial  is  fitted  to  the  data  points  by  any  standard  least  squares 
curve  fitting  technique. 

In  particular,  if  a  series  of  measurements,  X,,  X2,...,Xn  of  a  parameter  such  as  a  rectangular 
coordinate  are  made  as  a  function  of  time  and  a  polynomial  of  degree  p  is  fitted  to  this  data  and  if  the 
corresponding  points  on  the  curve  are  X',,  X2  then: 


i-1 


(4.5.1) 


is  an  estimate  of  the  standard  deviation  of  the  random  error  in  the  measurements  XI  ,X2...,Xn. 

If  p,  the  degree  of  the  polynomial,  is  small  compared  to  n,  the  number  of  points,  the  value  of  ox 
approximates  the  usual  definition  of  RMS  error,  i.e., 


RMS  - 


(x1  -  xp2 


(4.5.2) 


Estimated  standard  deviations  are  frequently  plotted  against  the  time  of  the  mid-point  of  the  span, 
and  significant  information  concerning  the  variation  of  random  errors  with  time  or  geometry  can  be 
obtained  from  these  plots. 


If  the  error  is  serially  correlated  then  the  estimates  of  the  variance  of  the  random  error  will  be 
biased  depending  on  the  extent  of  the  correlation.  This  is  discussed  in  more  detail  in  Section  4.7. 

However,  in  estimating  the  errors  by  either  the  variate  difference  and/or  least  square  curve  fitting 
methods,  it  is  usually  assumed  that  autocorrelation  and  cross  correlation  are  zero  or  very  nearly  so.  In 
actuality,  the  following  possible  conditions  may  exist  in  data  acquired  by  range  instrumentation: 

1.  The  errors  in  the  measurements  X,  Y,  Z  at  each  point  in  time  are  not  independent-cross 
correlation, 

2.  The  errors  in  successive  trajectory  points  are  not  mutually  independent-serial  or 
autocorrelation, 
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3.  All  measurements  of  a  given  type  do  not  possess  a  common  variance  heteroscedastic 
(versus  homoscedastic), 

4.  Time  measurements  are  not  error  free-time  error. 

Whenever  one  or  all  of  the  above  conditions  exist  in  the  acquired  data,  more  advanced  and 
sophisticated  techniques  are  needed  to  estimate  the  random  errors  in  the  data. 
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4.6  USE  OF  A  DIGITAL  HIGH  PASS  FILTER 


In  Section  3.0  random  error  was  characterized  in  the  frequency  domain  as  comprising  the  higher 
frequency  portion  of  the  power  spectrum  of  the  error  distribution.  If  this  is  the  case,  then 
frequency-constraining  digital  Filters  can  be  constructed  which  will  effectively  separate  the 
high-frequency  components  from  the  low  at  a  designated  cut-off  frequency.  This  technique  is  often 
more  desirable  than  using  a  least  square  polynomial  to  smooth  the  data  and  then  subtracting  the 
smooth  data  from  the  raw  measurement.  The  reason  for  the  desirability  of  the  digital  Filter  over  the 
least-square  polynomial  is  that  the  digital  filter  can  be  designed  so  that  its  frequency  response  has  a 
sharper  roll  off  with  less  “leakage"  or  side  lobe  effect  so  that  undesirable  frequencies  do  not  pass 
through  the  Filter.  Examples  of  second-degree  polynomial  frequency  responses  over  31  and  21  points 
and  digital  Filters  with  approximately  the  same  cut-off  frequencies  are  given  in  Figure  5.  It  is  noted 
from  the  graph  that  the  high-frequency  “leakage”  of  the  polynomial  fits  is  as  high  as  25  per  cent.  The 
corresponding  leakage  for  a  low-pass  digital  Filter  over  25  points  with  a  cutoff  of  0.05  cycles  per  unit 
time  is  2  per  cent. 

When  the  residual  distribution  is  obtained  from  a  high-pass  numerical  filter,  the  standard 
deviation  and  RMS  of  the  distribution  can  be  estimated  using  relations  (4.5.1)  and  (4.5.2)  where  p  is 
the  highest  degree  of  a  polynomial  which  will  pass  through  the  corresponding  low-pass  filter  without 
change. 
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EFFECTS  OF  AUTOCORRELATION  OF  THE  ESTIMATION  OF  RANDOM  ERRORS 


In  Section  3.3  it  was  pointed  out  that  high  autocorrelation  makes  the  random  error  appear 
smoother  and  can  lead  one  to  false  conclusions  concerning  reduction  techniques.  This  is  due  to  the 
fact  that  high  correlation  may  appear  as  signal  information.  Because  of  this  trait,  autocorrelated 
random  errors  will  bias  the  data,  and  the  amount  of  bias  usually  depends  on  how  high  the  correlation 
is. 


To  illustrate  how  autocorrelation  affects  data  reduction  and  random  error  estimation,  the 
different  random  errors  illustrated  in  Figure  2  were  superimposed  on  a  known  quadratic  function,  a0  + 
aj  t  +  a2 12 ,  with  a0  =  -1,  at  *  -1,  &2  =  1 »  and  the  variance  of  each  error  distribution  was  o2  =  0.01.  The 
results  of  approximating  the  quadratic  signal  with  uncorrelated  random  error  are  shown  in  Figure  6 
while  the  results  with  high  autocorrelated  error  are  given  in  Figure  7.  Each  of  these  figures  displays  a 
complete  picture  of  the  analysis.  The  First  graph  displays  the  predicted  curve  (line)  with  the  raw  data 
(signal  plus  noise)  on  the  same  graph  as  circles.  The  second  graph  displays  the  residuals.  The 
information  on  the  bottom  of  the  figures  describes  how  good  the  approximation  is.  A  comparision  of 
the  two  figures  is  summarized  in  the  following  table 


TABLE  I 


Auto 

Percent  Error 

Standard  Error 

Variance 

Best  Degree 

Correlation 

(rho) 

a0 

®.i 

*2 

c 

a0 

0 

0 

Estimate 

Fit 

0 

1.852 

3.808 

1.474 

.07311 

.04750 

.02310 

.0095 

2 

.99 

20.687 

29.893 

11.406 

.05548 

.02418 

.01176 

.0025 

10 

Table  I  shows  that  highly  correlated  error  gives  invalid  results  in  every  area.  The  error  in  the 
regression  coefficients  with  rho  =  .99  is  much  higher  and  the  estimates  of  the  standard  errors  in  the 
regression  coefficients  are  grossly  underestimated  and  are  not  consistent  with  the  actual  errors.  With 
rho  =  0  the  variance  estimate  is  very  close  (.0095  versus  0.01)  while  the  variance  estimate  for  rho  =  .99 
is  underestimated.  It  is  noted  that  the  best  polynomial  fit  was  correct  with  uncorrelated  noise. 
However,  when  rho  =  .99  was  present,  the  best  fit  was  a  polynomial  of  degree  10.  The  example  given 
has  shown  that  data  containing  random  error  with  high  autocorrelation  can  cause  the  following 
adverse  effects. 

1.  Often  makes  data  appear  smoother. 

2.  May  appear  as  signal  or  information. 

3.  Biases  estimates  of  regression  coefficients. 


4. 


Biases  estimates  of  variance  of  error  distribution. 


5.  Biases  estimates  of  variance  of  regression  coefficients. 

6.  Biases  “goodness-of-fit”  tests. 

The  effects  of  correlated  error  on  estimates  of  the  variance  when  the  variate  difference  method  is 
used  is  also  adverse.  High  correlation  causes  the  variance  estimates  to  be  low.  If  the  autocorrelation  is  a 
first  order  variate  difference  on  each  point  to  eliminate  a  quadratic  trend,  it  can  be  shown  that  the 
variance  estimate  will  be 


82  -  (1  -  p)(l  -  . 5p  +  .lp2)  a2 


(4.7.0) 


where 


p  -  correlation  coefficient  of  lag  1 
a 2  *  variance  of  error  distribution 

If  we  use  the  same  technique  on  every  other  point  (i.e.,  with  a  lag  of  2)  then  the  variance 
.stimate  is 


62  -  (l  -  o2)(l  -  .8p  +  .2a2)  a2 


(4.7.1) 


The  results  of  (4.7.0)  and  (4.7.1)  are  illustrated  in  Figure  8.  Note  from  the  figure  that  with  a  lag 
of  1  (every  point)  and  with  p  =  0.9,  the  variance  of  the  error  is  underestimated  by  a  factor  of  0.0631. 

It  is  stressed  that  there  are  techniques  which  may  be  used  to  estimate  and  compensate  for 
autocorrelated  error.  If  the  correlation  is  known  a  priori,  it  can  be  compensated  for  by  “extended”  or 
weighted  least-square  techniques  which  utilize  the  inverse  of  the  autocovariance  matrix  to  obtain  a 
weighted  estimate  which  compensates  for  correlation.  If  the  autocovariance/correlation  is  not  known, 
it  can  be  approximated  by  analysis  of  residuals  from  first  estimates.  An  iterative  approximation 
process  can  be  set  up  to  obtain  updated  estimates  of  the  autocovariance  functions  which  are  used  for 
updated  estimates  of  the  regression  coefficients.  Information  regarding  the  autocorrelation  function 
may  often  be  obtained  by  establishing  a  history  of  data  analysis  from  previous  operations. 
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FIGURE  7  Leost  squore  regression  on  y(t)  =  -  I  -  t  +  r  +  Nt 
Nf  =  .99  Nt.(  +  ijj  is  highly  correlated  with  p  =  .99,  mean 


=  0  and 


variance 


=  0.01 


FIGURE  8 


REDUCTION  FACTOR  K  IN  ESTIMATION  OF 
O’2  USING  3rd  ORDER  VARIATE  DIFFERENCE 
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5.0 


INTERVAL  ESTIMATION  OF  RANDOM  ERRORS 


In  the  previous  section,  techniques  for  point  estimation  of  the  variance,  standard  deviation  or 
root  mean  square  of  the  random  error  distribution  were  discussed.  It  should  be  noted  that  the 
estimates  of  these  parameters  are  themselves  random  variables  and  have  corresponding  probability 
distribution  functions.  This  fact  points  out  that  the  use  of  sample  random  variables  (or  residuals  in  the 
case  of  error  analysis)  involves  the  concept  of  confidence,  tolerance  and  prediction  regions.  In  the  case 
of  one  dimensional  random  variables,  these  regions  are  intervals  on  the  real  line.  The  regions 
corresponding  to  higher-dimension  random  variables  will  be  discussed  in  Section  8.0  after  the 
covariance  matrix  has  been  discussed. 

Interval  estimation  is  an  important  concept  and,  in  general,  if  we  are  given  a  sample  of  n  residuals 
ei...,en  from  the  population  e,,...,en,  the  problem  is  to  find  two  statistics  (L,  U)  defining  an  interval 
which  has  100a%  probability  of  containing  100t%  of  the  individual  values  in  the  population.  Such  an 
interval  is  called  a  tolerance  interval  with  lower  and  upper  tolerance  limits  L  and  U.  The  value  of  a,  0 
<  a  <  1,  is  called  the  confidence  coefficient.  For  example,  if  a  =  .90  and  y  =  .95,  we  are  90% 
confident  that  at  least  95%  of  the  individual  values  will  be  in  the  interval  (L,U). 

Besides  interval  estimation  for  individual  values  of  the  population,  we  are  interested  in  interval 
estimates  for  the  parameters,  0,  which  define  the  probability  density  function  of  the  population.  We 
want  to  find  statistics  L(0)  such  that 

Prob  (L(0)  <  0  <  U(0))  =  y  (5.0.1) 

where  0  would  be  a  parameter,  such  as  the  mean  or  variance  of  the  distribution.  In  this  case, 
the  interval  L(0),  U(0)’is  called  a  confidence  interval  for  the  true  (but  unknown)  value  of  the  parameter. 
Another  interval  estimate  is  the  prediction  interval.  In  this  case,  we  construct  an  interval  that  has  the 
preassigned  probability  of  containing  the  next  observation. 


Preceding  page  blank 
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5.1  THE  CASE  OF  THE  NORMAL  DISTRIBUTION 


It  is  often  assumed  that  the  random  error  is  a  random  variable  which  has  the  normal  or 
Gaussian  distribution.  This  assumption  may  be  based  on  experience  or  used  because  it  is  well  known, 
convenient,  and  gives  approximately  correct  results.  It  is  pointed  out  that  in  the  real  world  many 
random  variables  are  not  normally  distributed,  and  this  includes  random  errors.  An  example  of  an 
error  which  is  uniformly  distributed  (rather  than  normally  distributed)  is  round-off  error.  However, 
the  importance  of  the  normal  distribution  in  error  analysis  is  recognized  and  some  of  its  basic  uses  arc 
considered  in  this  document. 

A  Figure  of  the  standard  normal  distribution  (i.e.,  M  =  0,  a2  =  1)  is  given  in  Figure  9.  If  ct  is  a 
sample  from  the  population  et  which  corresponds  to  the  normal  distribution  with  mean  n  and  variance 
o2  (assumed  known  in  this  instance),  then  tolerance  intervals  of  the  form  fit  ka  corresponding  to  a  = 
1  (absolute  certainty)  can  be  constructed.  The  values  of  k  corresponding  to  the  more  widely  used 
values  of  7  are  given  in  the  following  table 

AREAS  (1007%)  UNDER  THE  NORMAL  CURVE  IN  THE  INTERVAL  n±ka 
k  7 


.6745 

.5000 

.7979 

.5751 

1.0000 

.6827 

1.6449 

.9000 

2.0000 

.9545 

3.0000 

.9973 

These  particular  values  of  k  and  7  correspond  to  the  following  error  quantities 

.6745o:  probable  error  (PE) 

.7979 0:  mean  absolute  error  (ME) 
a:  standard  deviation 
1.6449a:  map  accuracy  standard  (MAS) 

2o:  significant  error 

3a:  highly  significant  or  near  certainty  error 

Another  important  quantity  is  the  standard  error  of  the  mean,  or  standard  error  of  estimate.  If 
the  mean  and  variance  of  a  sample  are  estimated  from  n  points,  then  he  standard  error  is  a/J n. 
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RELATIVE  FREQUENCY 


Y  -  —  e_X2/2,  u  -  0,  0  -  1 
X  -  Z/o 

a  -  Standard  Deviation 
Z  -  Actual  Deviation 


FIGURE  9.  Standard  Normal  Frequency  Distribution  Curve. 
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5.2  CONFIDENCE  INTERVALS  ON  ERROR  PARAMETER  ESTIMATES 


Suppose  ej,...,en  is  a  sample  observed  from  the  population  et  ~  normal  (jU,o2)  where  the 
population  parameters  n  and  a2  are  unknown.  Assuming  that  the  observations  are  independent,  then  /u 
and  a  are  estimated  by  the  sample  mean  F  and  sample  standard  deviation  se,  where 


n 


i-1 


(5.2.1) 


■o  ■  J ^rr  l  <*1  -  •> 

’  i-1 


F  and  se  are  random  variables  with  e  ~  normal  (n,a2  /n)  and  the  quantity 


(n  -  1)«* 

Z - 1 


(5.2.2) 


has  the  Chi-square  distribution  with  n  -  1  degrees  of  freedom  (i.e.,  Z  has  X2  ,  distribution).  It  can 
further  be  shown  that  the  statistic 


t  .  £  ~  u  (5.2.3) 

has  the  “student”  T-distribution  with  n  -  1  degrees  of  freedom.  From  relation  (5.2.3)  we  can  obtain  a 
confidence  interval  for  n  by 


Prob 


t 


l  n-‘ 


j 


-  Y 


or 


e  ±  t 

n-i 


(5.2.4) 
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is  a  100t%  confidence  interval  on  n,  where  t  depends  on  the  confidence  coefficient  a  and  the  degrees 
of  freedom  v  which  in  the  case  of  n  observations  correspond  to  (n  -  1).  When  tvot  =  1,  then  the 
confidence  interval  is  "e  ±se//rT  and  the  quantity  se/Jn  is  called  standard  error  of  the  mean  or  just 
“standard  error.”  Values  for  ty  a  can  be  obtained  from  tables  in  most  statistical  texts. 

In  a  similar  manner  we  can  get  a  confidence  interval  on  0  utilizing  (5.2.2).  From  (5.2.2)  we  get 


( 


Prob 


< 


<n  -  U  sj 


-  Y 


(5.2.5) 


which  leads  to  a  100 7%  confidence  interval  on  a2  which  is 


(n  -  1)b2  (n  -  l)s2 

- y  —  —  <  O2  <  7~~*~ 

X(^)  X(^) 


(5.2.6) 


Values  of  can  be  obtained  from  Chi-square  tables  in  statistical  texts. 
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6.0  METHODS  USED  TO  EVALUATE  SYSTEMATIC  ERRORS 


The  problem  of  evaluating  systematic  errors  in  a  given  instrumentation  system  is  considerably 
more  difficult  than  that  of  evaluating  random  errors,  unless  a  suitable  comparison  standard  is  available. 
Unfortunately,  the  systematic  errors  can  be  several  orders  of  magnitude  greater  than  the  random  errors 
and  no  absolute  comparison  standard  exists.  There  are,  however,  several  approaches  which  can  be  used 
to  evaluate  systematic  errors. 

6.0.1  ANALYTICAL  INVESTIGATION  OF  POSSIBLE  ERROR  SOURCES 

This  involves  representing  the  system  by  as  comprehensive  a  mathematical  model  as  possible. 
Reasonable  perturbation  of  each  of  the  model  parameters  is  introduced  and  the  resultant  effects  of 
this  on  the  data  are  calculated.  From  this  analysis  may  be  obtained  an  estimate  of  the  likely  range  of 
systematic  errors.  This  approach  is  limited  by  the  store  of  scientific  knowledge,  which  makes  it 
extremely  difficult  to  construct  a  sufficiently  comprehensive  mathematical  model  of  the  system.  This 
would  appear  to  be  especially  true  of  propagation  anomalies  which  particularly  affect  electronic 
systems. 


6.0.2  COMPARISON  OF  OBSERVATIONS  FROM  DIFFERENT  SYSTEMS 

If  it  is  known  that  System  I  is  significantly  more  accurate  in  an  absolute  sense  than  System  II, 
then  an  estimate  of  the  systematic  error  in  II  can  be  made  by  noting  the  discrepancies  between 
observations  of  the  two  systems.  The  estimate  of  the  total  error  in  System  II  is  made  by  computing 
tl  e  difference 


A  =  (System  II  -  System  I) 

at  each  time  point.  The  difficulty  with  this  approach  lies  in  selecting  a  suitable  standard.  Emphasis 
must  be  placed  upon  the  fact  that  upper  bounds  on  possible  systematic  errors  in  the  standard  must  be 
known,  and  these  bounds  must  be  narrow  compared  with  reasonably  lower  bounds  for  the  systematic 
error  in  the  system  being  investigated.  Thus,  for  the  calibration  standard  it  is  absolute  accuracy,  not 
relative  or  internal  accuracy,  which  is  of  prime  importance.  This  fact  makes  it  unprofitable  to  compare 
systems  indiscriminately.  Also,  it  is  necessary  to  employ  a  mathematical  model  of  the  system  being 
treated  in  order  to  transform  the  standard  data  and  propagate  any  error  in  the  standard  into  the  tested 
system. 

The  observations  collected  and  used  for  error  analysis  may  be  from  actual  operational  tests  or 
from  tests  which  are  specifically  designed  for  the  purpose  of  error  analysis.  Tests  of  the  latter  type  are 
called  calibration  tests.  The  most  simple,  economical,  and  frequently  used  calibration  technique  is  the 
“static”  'calibration  test.  This  utilizes  calibration  equipment  and  objects  as  “line-ups,”  boresight 
towers,  and  stars.  In  most  instances,  however,  the  static  calibration  test  is  not  effective  in  estimating 
the  many  errors  encountered  during  an  actual  operation.  The  reason  for  this  is  that  most  of  the 
systematic  error  components  are  dynamic  in  nature  and  depend  upon  variables  such  as  distance, 
tracking  rates,  atmospheric  conditions  (temperature  pressure,  humidity,  etc.),  and  the  direction  in 
which  the  instrumentation  system  is  pointing.  As  a  result,  dynamic  calibration  tests  utilizing  a  good 
test  design  and  a  standard  with  sufficiently  small  error  provides  a  better  estimate  of  the  systematic 
error. 
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6.0.3  ANALYSIS  OF  MEASURING  RESIDUALS  FROM  OVER-DETERMINED  LEAST 
SQUARES  SOLUTIONS 


When  a  given  instrumentation  system  provides  a  redundancy  of  data,  a  least  square  solution  may 
be  performed.  The  adjustment  leads  to  a  set  of  measuring  residuals  (or  corrections)  for  each  trajectory 
point.  The  criterion  for  adjustment  is  the  principle  of  maximum  likelihood.  If  these  residuals  indicate, 
from  point-to-point,  a  randomness  about  zero,  there  is  good  indication  that  the  systematic  error  is 
small.  If  the  residuals  do  not  display  this  randomness  about  zero,  one  must  conclude  that  the  original 
observations  are  biased  or  have  systematic  error.  A  total  error  computed  from  the  residuals  gives  an 
estimate  of  the  upper  limit  of  the  systematic  error  in  the  observations. 

6.0.4  ANALYSIS  OF  RESIDUALS  FROM  A  BEST  ESTIMATE  OF  TRAJECTORY 

When  there  is  simultaneous  track  with  more  than  one  instrumentation  system,  a  Best  Estimate  of 
Trajectory  (BET)  can  be  computed.  This  trajectory  and  the  residuals  of  the  individual  systems 
referenced  to  it  can  then  be  analyzed  to  determine  estimates  of  the  systematic  errors.  The  BET  is  a 
powerful  method  of  estimating  systematic  errors  in  both  the  reduced  data  and  the  measured 
parameters.  There  are  several  methods  of  obtaining  a  BET,  depending  upon  the  criteria  for  test.  The 
most  powerful  of  the  methods  developed  to  date  involves  the  use  of  mathematical  models  of  the  errors 
of  the  individual  systems.  In  all  the  applications  of  BET  to  error  studies,  the  basic  idea  is  that, 
according  to  the  criterion  agreed  upon,  a  best  estimation  of  what  the  observations  should  have  been  is 
made  and  then  compared  with  the  observed  data. 

6.0.5  MULTILATERATION 

As  in  the  case  of  BET,  when  redundant  data  from  several  instruments  is  available  a 
multilateration  technique  which  uses  inputs  from  several  different  sources  such  as  doppler,  pulse  radar 
data,  and  inertial  guidance  data  can  be  used.  This  technique  has  the  error  model  design  as  part  of  the 
solution  and  attempts  to  “weight  out”  the  effects  of  error  sources  by  associating  each  source  with  an 
uncertainty  which  is  estimated  from  a  priori  testing  and  is  updated  by  the  multilateration  solution.  If 
error  sources  are  correlated,  they  adversely  affect  the  solution  and  must  be  taken  out,  using  estimates 
from  calibration  performed  prior  to  the  test.  Once  this  is  done,  estimates  of  the  magnitudes  of  all  the 
error  model  coefficients  are  made  simultaneously  employing  a  maximum  likelihood  technique,  which 
is  mechanized  using  the  Kalman  filter  approach. 


J 
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6.1 


ERROR  MODELS 


In  Section  3.0  it  was  stated  that  systematic  errors  are  deterministic.  This  means  the  systematic 
error  S(t)  can  be  represented  analytically  in  terms  of  a  mathematical  model.  Since  the  random  error  is 
nondeterministic,  it  must  be  represented  by  some  probability  model.  In  general,  one  may  set  up  error 
models  which  fall  into  the  following  categories: 

I.  Deterministic 

Parametric  estimation— an  analytic  model  which  describes  the  error  in  terms  of  parameters 
which  are  the  components  of  the  error. 

Nonparametric  estimation— a  numerical  model  which  approximates  the  error  numerically, 
i.e.,  as  a  sequence  of  numbers  such  as  the  output  from  a  digital  filter. 

II.  Nondeterministic 

Parametric  estimation— a  stochastic  equation  depending  upon  parameters  which  are  random 
variables,  e.g.,  the  expression  of  a  random  process  by  a  pth  order  Markov  Process.  If  the  process  is 
stationary  it  can  be  characterized  in  terms  of  parameters  which  index  its  probability  density  function. 

The  most  simple  systematic  error  model  would  be 

S(t)  =  /i  (6.1.1) 

In  this  case  the  systematic  error  is  a  constant  bias,  “offset”  or  “zero  set”  error.  If  S(t)  is 
estimated  by  a  sequence  of  n  Aj’s,  then 

n 

y  *  S  -  i  l  Aj  (6.1.2) 

j-1 

The  Aj’s  also  contain  random  error  and  the  averaging  in  (6.1.2)  separates  ju  from  N(t).  The 
estimate  of  the  random  error  parameter  would  be 

9  -  14  W-®*.  <6-1'31 

j-1 

If  the  error  components  are  independent,  then  the  error  model  can  be  expressed  as  a  linear 
combination  of  independent  variables,  i.e., 

k 

Gi  "  l  °jXi  +  ni  (6.1.4) 

j-1 
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The  error  model  given  in  (6.1.1)  can  be  expanded  to  include  terms  which  describe  various  error 
components.  A  good  example  is  a  typical  error  model  frequently  used  for  the  AN/FPS-16  radars.  The 
error  models  for  range,  azimuth,  and  elevation  are  given  below. 


AR 


•j  ■  ao  +  aiRj  +  a2Rj  +  asRj  +  ai,f(S/N)  +  asRj  +  a6RjRj  (6.1.5) 


+  a7  esc  Ej  +  a8  sin  R.  +  a9  cos  7777  R,  +  a 


6000  "j 


y i  zi 

+  *ll  p  ai2  p  +  N, 

Rj  J 


6000  j 


10 


21 


where 


a0  =  constant  bias 
at  =  timing  delay 
a2  =  acceleration  servo  lag 
a3  =  jerk  servo  lag 
a4  =  beacon  delay 
a$  =  oscillator  drift  or 
scale  factor 
a6  =  time  dilation 


a 7  =  residual  refraction 
a8  =  resolver  nonlinearity 
a9  =  resolver  nonlinearity 
a,  o  =  x  survey  error 
a  i ,  =  y  survey  error 
a  j  2  -  z  survey  error 

Nj  *  random  error 


AAj  -  b0  +  b^  +  b2Aj  +  b3Aj  +  b9sin  Aj  tan  Ej  +  b5  cos  A^  tan  E^ 
+  bsRjAj  +  b7  sec  Ej  +  b a  tan  E^  +  b9  sin  Aj  +  b10  cos  Aj 
+  ba  —  cos  A.  sec  E  +  B12  sin  A.  sec  Ej  +  Nj 


(6.1.6) 


where 


b0  =  constant  bias 

b|  =  timing  bias 

b2  =  acceleration  servo  lag 

b3  =  jerk  servo  lag 

b4  =  mislevel 

bs  =  mislevel 

b6  =  time  dilation 


b7  *  collimation 
b8  =  nonorthogonality 
b9  =  encoder  nonlinearity 
bj  o  =  encoder  nonlinearity 
b]  i  =  x  survey 
b|  j  =  y  survey 
Nj  =  random  error 


(6.1.7) 


AEj  ■  co  +  ciEj  +  C2Ej  +  cjEj  +  cn  tin  Aj  +  c$  cos  +  c 


+  07  cot  Ej  +  c o  sin  Ej  +  cj  cos  Ej  +  Cj  i  —  sin  Aj  sin  Ej 


+  ci 2  cos  A4  sin  Ej  +  ci 3  At  cos  Ej  +  N., 


\)  “j 


R(t) 


J  J 


where 


c0  =  constant  bias 

c,  =  timing  delay 

c2  =  acceleration  servo  lag 

c3  -  jerk  servo  lag 

c4  =  mislevel 

cs  =  mislevel 

c6  =  time  dilation 


c7  =  residual  refraction 
c8  =  encoder  nonlinearity 
c9  =  encoder  nonlinearity 
c,  o  =  x  survey 
Cj  1  =  y  survey 
C|  2  =  z  survey 
Nj  =  random  error 


In  the  linear  error  models  the  dependent  variable  is  the  total  error  residual,  and  the  independent 
variables  are  usually  assumed  to  be  uncorrelated  and  are  determined  by  the  physics  or  geometry  of  the 
situation.  The  independent  variables  are  assumed  to  be  known  without  (or  at  least  with  negligible) 
error.  Since  the  independent  variables  depend  on  the  physics  and  geometry  of  a  particular  test  for  a 
certain  instrumentation  system,  it  is  important  to  stress  that  the  error  model  depends  not  only  on 
analysis  of  the  physics  involved,  but  on  the  design  of  the  test  itself.  For  example,  if  one  had  a 
near-perfect  standard  for  comparison  for  obtaining  good  error  estimates,  one  could  not  obtain  a  valid 
estimate  of  the  mislevel  error  coefficients  if  the  test  design  failed  to  enable  the  instrumentation  system 
to  traverse  less  than  180°  in  azimuth.  Also,  one  could  not  estimate  scale  factor  error  if  the  range 
remained  approximately  constant.  It  is  obvious  that  in  order  to  estimate  dymanic  errors  one  must  have 
a  corresponding  dynamic  test  design. 

One  proposed  test  for  dynamic  error  estimation  involves  the  use  of  a  calibration  satellite.  This 
type  of  vehicle  would  allow  tracking  in  virtually  every  quadrant  over  a  wide  variety  of  rates,  look 
angles,  and  positions.  Further,  its  trajectory  can  be  predicted  very  accurately  using  the  equations  of 
motion  as  a  model.  The  orbital  parameters  can  be  determined  by  BET  techniques,  using  world-wide 
networks  of  instrumentation  systems.  Many  problems  would  still  exist  in  a  program  of  this  type,  such 
as  unmodeled  parameters  in  the  orbital  model,  nonuniformity  in  the  earth’s  gravitational  field,  etc. 
However,  during  certain  satellite  passes  and  combinations  of  passes,  the  effects  of  such  errors  can  be 
minimized.  For  further  detailed  information  on  the  use  of  a  calibration  satellite,  see  Reference  24. 
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6.2  REGRESSION  ANALYSIS 


The  linear  error  models  discussed  in  the  previous  section  are  generally  approximated  utilizing 
multiple  linear  regression  techniques.  The  techniques  for  regression  analysis  are  discussed  in  this 
section  and  the  techniques  for  evaluating  the  validity  of  the  analysis  are  discussed  in  Section  6.3, 
Analysis  of  Variance. 

When  an  acceptable  standard  is  available  and  e,,  e2,  ...,  en  obtained  they  are  assumed  to  be  a 
realization  from  the  process 

ci  ■  “o  +  “‘Xll  +  “2*i2  +  •"  +  “kXlk  +  nl  l6'2’1) 


i  =  1,  2,  ...,  n.  The  a’s  are  unknown  parameters  and  is  the  random  error,  tj  is  assumed  to  be  a 
stochastic  process  whose  elements  are  homoscedastic,  independent,  with  mean  =  o  and  variance  o*  . 
In  matrix  form  (6.2.1)  is 


e-Xa  +  rj 


(6.2.2) 


where  e  and  are  n  x  1  column  vectors,  a  is  a  (k-*  1)  x  1  column  vector  and  X  is  an  n  x  (k+1)  matrix  of 
known  elements 


1 

*11 

X12 

. . .  Xlk 

1 

• 

x21 

• 

x22 

• 

. . .  X2|[ 

•  • 

• 

• 

1 

• 

• 

Xnl 

• 

• 

X*l2 

•  • 

•  • 

Xnk 

(6.2.3) 


The  matrix  X  consists  of  the  n  values  of  the  k  independent  variables  and  is  called  the  design 
matrix.  In  the  analysis  and  design  of  experiments,  the  independent  variables  are  called  factors,  the 
dependent  variable  (in  this  case  the  error)  is  called  the  response  variable  or  yield.  The  n  values  of  the 
independent  variables  are  called  levels. 


In  the  standard  least  squares  analysis  of  (6.2.2),  estimates  a  are  made  of  the  parameters  a  such 
that  the  sum 


n 

S  ^  "  *o  *lXil  ” 

i-1 


Vik>2 


is  minimized.  This  done  by  taking 


(6.2.4) 


3S  _  3S  _  .2£.o. 

3*0  3«1  •"  *k 

Preceding  page  blank 


(6.2.5) 
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(6.2.6) 


The  system  (6.2.5)  gives  the  normal  equations  for  a.  These  are 
(XTX)a  =  XTe . 

and  it  can  be  shown  that  (X^X)  is  a  positive  definite  matrix  and  has  an  inverse.  Thus 

a=(XTXpXTe  (6.2.7) 

The  a’s  are  called  regression  coefficients  and  are  least  square  estimators  of  the  parameters  a.  Since 
the  a’s  are  dependent  upon  the  ej’s  they  themselves  are  random  variables  and  are  characterized  by  a 
kth  dimensional  probability  density  function.  Since  the  system  of  normal  equations  was  linear,  the 
variance  of  the  regression  coefficients  is  estimated  by  the  sample  covariance  matrix: 


where 


Var(a)  -  S2  -  [XTX]_1s2 


n  N  2 

I  rdhr-p 

i-l 


(6.2.8) 


(6.2.9) 


Sj(j  is  the  estimate  of  the  variance  of  rj,  computed  from  the  residuals  Nj,  where 


N 


(6.2.10) 


i.e.,  the  N’s  are  the  residuals  about  the  regression  estimate  of  e. 


If  the  tj’s  are  assumed  to  be  characterized  by  the  normal  distribution,  then  since  the  a’s  are  linear 
combinations,  they  are  also  normal  variates  with  expected  value  estimated  by  (6.2.7)  and  variances 
estimated  by  (6.2.8).  Confidence  intervals  on  the  parameters  o^.  and  cq  can  be  obtained  using  the 
Chi-square  and  normal  distributions  as  discussed  in  Section  5.2. 


If  the  errors  are  not  independent,  then  it  can  be  shown  that  if  this  is  not  accounted  for,  the 
estimate  of  a  will  be  biased.  The  techniques  for  handling  correlated  error  come  under  regression 
analysis  with  “extended”  or  weighted  least  squares  approximation.  A  detailed  discussion  of  this  topic 
is  found  in  Reference  26. 
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6.3  ANALYSIS  OF  VARIANCE 


Many  of  the  previous  sections  show  that  error  analysis  is  actually  the  analysis  of  the  variation  of 
measurements  about  some  accepted  standard.  The  problem  involved  in  error  analysis  is  to  evaluate  it 
in  terms  of  explained  variation  (systematic  error)  and  unexplained  variation  (random  error).  If  one 
should  obtain  error  estimates  from  several  sources  such  as  various  tests,  different  instrumentation 
sites/systems,  several  runs,  flights,  times,  aircraft,  missiles,  etc.,  he  might  desire  to  consolidate  the 
results  and  design  a  technique  to  break  the  total  variation  into  components  which  will  explain  the 
effects  these  varied  sources  have  on  the  total  "ariation.  For  example,  the  error  models  in  Section  6.1 
give  several  sources  which  will  contribute  t  .  le  total  variance  of  the  error.  The  questions  are  “how 
much  do  they  contribute”  and  “is  the  contribution  due  to  any  one  source  significant?”  If  the  answers 
to  questions  such  as  these  can  be  obtained,  we  would  then  know  which  variables  to  use  or  delete  in 
our  regression  analysis.  The  technique  used  to  break  variation  in  data  into  source  components  is  called 
Analysis  of  Variance,  which  is  highly  dependent  upon  the  design  of  experiments  briefly  mentioned  in 
Section  6.2. 

The  sources  of  variation  considered  in  analysis  of  variance  are  called  variables  or  factors.  The 
factors  may  be  quantitative  (such  as  a  dynamic  error  model)  or  qualitative  (such  as  different  missiles 
or  radars).  Analysis  of  variance  can  correspond  to  several  test  designs.  The  most  simple  test  design  is 
the  “completely  randomized”  or  between-group  and  within-group  design.  In  this  case  there  is  only  one 
source  of  variation.  For  example,  suppose  that  total  error  estimates  y  are  made  for  a  particular  radar 
for  k  different  tests  from  the  same  type  of  missile.  The  data  would  be  grouped  according  to  tests  and 
the  problem  is  to  estimate  the  amount  of  variation  (out  of  the  total  variation)  between  tests;  the 
amount  of  variation  within  tests,  and  to  test  whether  or  not  such  variation  is  “significant.”  Such  a  test 
is  called  completely  randomized  because  the  grouping  of  the  tests  is  accomplished  in  a  random  manner 
to  eliminate  any  systematic  trends  which  might  exist  from  test  to  test.  The  mathematical  model  for 
such  a  design  is 


' "  +  “i  +  Eir  3  ‘  1 . V  1 ' 1 . k  (6'3,1) 

where 

yy  =  ith  measurement  from  the  jth  group 

nj  =  number  of  observations  in  the  ith  group 

li  =  mean  of  population  from  which  groups  are  sampled 

otj  =  deviation  of  mean  of  ith  group  from  n 

Cjj  =  unexplained  variation  in  yy. 
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The  a’s  are  regarded  as  random  variables  with  mean  zero  and  variance  .  The  term  Cj:  is  often  called 
experimental  error,  and  in  the  between-within  experiment  o|j j  would  be  the  variance  within  the  groups. 
Table  6.3.0  gives  a  typical  between-within  analysis  of  variance  design  (one-way  analysis  of  variance) 
and  illustrates  how  the  parameters  are  estimated. 

TABLE  6.3.0 


ONE  WAY  ANALYSIS  OF  VARIANCE 


If  the  errors  are  not  independent,  the  analysis  of  variance  is  not  valid.  Errors  from  different 
groups  may  reasonably  be  assumed  uncorrelated,  but  errors  within  groups  are  often  correlated, 
especially  if  the  measurements  within  a  group  come  from  a  time  series.  However,  this  can  be  offset  by 
randomizing  within  each  group.  Thus,  in  order  to  be  effective,  the  one-way  analysis  of  variance  should 
be  completelyrandomized. 


The  analysis  of  variance  usually  has  a  “test  of  significance”  associated  with  it.  That  is,  we 
hypothesize  that  the  means  of  each  group  are  equal,  then  (under  the  assumption  the  groups  have  been 
completely  randomized  and  are  independent  normal  variates),  test  for  the  rejection  of  this  hypothesis 
at  a  prescribed  probability  level.  The  most  widely  used  test  is  the  ratio  of  the  variance  from  the  source 
(between)  to  the  variance  of  the  experimental  error.  This  forms  a  statistic  which  has  the  F  probability 
distribution  with  (k-1)  and  (N-k)  degrees  of  freedom.  This  is  the  same  as  saying 

K 

FfN^jdF.  (6.3.2) 

Jo  (k-i) 

Suppose  we  select  K  such  that  Prob(Sf /S^  <  k)  =  .95,  and  then  set  up  the  hypothesis  H0:  Mi  =  Mt 
=’"~  ^  means  each  group  are  equal).  From  the  analysis  of  variance  we  estimate  F 

-  Sf 

F  ■  sf  •  (6.3.3) 


1  ^  variance  between  groups  <  ^ 
Sf  variance  within  groups 
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A. 


if  F  is  then  compared  to  the  theoretical  value  of  F  from  the  F  distribution,  and  F  >  F  we  reject  H0  and 
say  the  group  means  are  significantly  different.  If  ¥  <  F  we  cannot  reject  H0  and  there  is  no 
significant  difference. 

When  certain  factors  which  are  known  are  placed  in  the  experimental  design,  additional 
constraints  are  added  to  the  design.  For  example,  we  could  select  the  p  radars  and  k  groups  from  the 
same  type  of  test.  We  block  the  k  groups  and  block  the  p  radars  in  each  of  the  k  groups.  This  type  is  a 
randomized  block  design  as  shown  below 
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The  analysis  for  the  randomized  block  design  is  called  a  two-way  analysis  of  variance,  and  the 
component  variances  measure  the  effects  between  tests  and  between  radars  with  respect  to  the 
variance  of  the  experimental  error.  If  the  cells  in  the  design  of  Figure  6.3.0  contain  r  replicates  each, 
then  it  becomes  possible  to  estimate  the  effects  of  interaction  between  radars  and  tests. 


The  completely  randomized  design  and  the  randomized  block  design  are  special  cases  of  a  more 
generalized  design  technique  called  factorial  design.  This  technique  handles  several  variables,  each  with 
several  levels  and  each  having  replicates.  This  type  of  analysis  could  be  used  in  conjunction  with  the 
error  models  (6.1.5)  and  (6.1.7).  The  error  model  for  radar  range  error  (Relation  6.1.5)  would  consist 
of  12  factors,  each  having  n  levels  (corresponding  to  n  observations),  with  one  replicate.  An  example 
of  an  analysis  of  variance  on  a  similar  radar  range  error  model  is  given  in  Table  6.3.1. 


Table  6.3.1  is  an  example  of  an  analysis  which  summarizes  Sections  6.0,  6.1,  6.2,  and  6.3,  and 
needs  some  explanation.  The  coefficient  of  determination  estimates  the  percentage  of  variation  out  of 
the  total  variation,  which  is  explained  by  the  regression  analysis.  The  computed  T  value  is  the  ratio  of 
the  regression  coefficient  to  its  standard  error  estimate.  The  partial  correlation  coefficient  is  an 
estimate  of  correlation  between  the  two  variables,  keeping  the  effects  of  the  other  variables  fixed.  The 
multiple  correlation  coefficient  is  the  correlation  between  the  observed  dependent  variable  and  the 
regression  estimate.  The  “proportion  of  variance  cumulated”  estimates  the  percent  of  the  total 
variation  due  to  each  individual  variable.  Examination  of  the  table  shows  that  the  fit  explains  only 
49%  of  the  total  variations  with  31.4%  and  17.2%  of  this  due  to  scale  factor  and  timing.  The  T  value 


shows  these  two  regression  coefficients  are  significantly  different  from  zero  at  the  95%  confidence 
level.  The  total  error  estimate  is  RMS  (AR)  =  9.7  feet.  The  bias  estimate  is  -7.54  feet,  while  the  noise 
estimate  is  4.63  feet.  Since  the  ratio  [(mean  AR)/d^)  ]  is  small  it  is  difficult  to  obtain  a  large 
coefficient  of  determination  (i.e.,  a  large  proportion  of  the  total  error  is  random  or  “unexplained”). 
However,  the  test  for  significance  and  T  values  indicates  the  regression  fit  is  significant. 
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TABLE  6.3.1 

ANALYSIS  OF  VARIANCE  OF  RADAR  ERROR  MODEL  (EXAMPLE) 
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53988  0.27293  -1.54071  1.09114  -1.41202  -0.08707  42.80281  0.00389 


7.0  THEORY  OF  ERROR  PROPAGATION 


In  addition  to  estimating  the  error  in  the  measurements  from  instrumentation  systems,  it  is 
necessary  to  estimate  the  error  in  the  final  reduced  data.  This  data  usually  consists  of  position, 
velocity,  and  acceleration  information  in  some  coordinate  (usually  rectangular)  system.  Thus,  errors  in 
range,  azimuth,  and  elevation  or  direction  cosines,  film  coordinates,  range  sums  or  differences  must  be 
propagated  into  rectangular  coordinates.  Furthermore,  if  velocity  and  acceleration  are  to  be  inferred 
from  this  positional  data,  the  manner  in  which  the  errors  propagate  must  be  known.  At  least  three 
factors  are  of  importance  here:  (1)  the  manner  in  which  observing  sites  are  located  with  respect  to 
each  other  (assuming  that  several  are  involved),  (2)  the  errors  in  the  measurement  systems,  and  (3)  the 
position  of  the  object  in  space  whose  coordinates  are  to  be  determined. 


Let  x  be  a  random  variable  with  mean  y  and  variance  a2.  If  y  ■  a0  +  *ix, 
it  can  be  shown  that 


g2  m  g2g2  (7 .0.1) 

y  1  * 

Further,  if  X} ,  X2,  X3,  . ...  xQ  are  independent  random  variables  with 
variance  a2,  then  the  sum 


n 

l 


y3  ’  V 


J  *  1 *2  » •  •  • 


(7.0.2) 


has  mean 


l 

i-l  Ji 


and  variance 


a|  o2. 

h  * 


(7.0.3) 


(7.0.4) 


Preceding  page  blank 
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N  is 


Further,  the  covariance  of  y  ^ ,  y^»  j,  k  ■  1,2 . 


f  ? 

■  /  a.  a . 

'jyk  h  ki 


(7.0.5) 


If  y  -  f(x)  is  some  function,  we  can  approximate  f(x)  in  the  neighborhood 
of  x  *  x  by  a  first  order  Taylor  series  expansion.  That  is 


y  ’  f(*o>  +  to  <*  -  V 
0 


(7.0.6) 


In  the  same  manner,  using  (7.0.1)  we  approximate  the  variance  in  y  by  the 
relation 


2  *  df  2  ? 
oz  ■  - —  a*  . 
y  dx  x 

1  o' 


(7.0.7) 


Now  suppose  y  «  f(xi,X2),  then  the  first  order  Taylor  Series  expansion  in  the 
neighborhood  of  (xi,X2)  •  (x°,X2)  is 


y  ;  f  M.*5)  +  ^  0  (*i  -  *1)  + |~  „  (<2  -  x2) , 


(7.0.8) 


and  from  (7.0.4)  and  (7.0.5)  the  form  of  o2  would  be 

2  faf  l2  2  j.  fiL.12  ,  ,  UL  „ 

°y  "  [3xJ  *1  13X2J  °x2  3xl  3x2  x1x2 


(7.0.9) 


general  if  y^  y^(xi»*2 . xn) »  ^  ■  l,2,...,m,  then  the  Taylor  Series 


expansion  estimate  of  the  covariance  matrix  of  y  is 


S2  -  RS  RT 

y  x 


(7.0.10) 
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(7.0.11) 


(7.0.12) 


(7.0.13) 


and 


nn 

R1  is  the  transpose  of  R 


If  the  xj’s  are  independent  random  variables  then  Sx  is  a  diagonal  matrix,  i.e.,  the  covariances  are 
all  zero.  In  this  case 
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i3xJ 

(7.0.14) 
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2  T 

The  relation  °y  "  *s  usuaUy  called  the  covariance  equation  and  is  classically  referred 

to  as  the  Gaussian  law  of  error  propagation.  This  technique  is  based  directly  upon  approximating  the 
change  in  the  dependent  variable  by  the  differential,  i.e.,  the  linear  portion  of  a  Taylor  Series 
expansion.  This  leads  one  to  wonder  how  much  error  is  involved  in  the  variance  estimates.  Actual 
examples  using  simulation  techniques  have  shown  the  method  to  be  very  good  and,  as  John  Tukey  is 
quoted  as  saying,  "The  most  important  conclusion  is  that  the  classical  propagation  formula  is  much 
better  than  seems  to  be  realized.  Examples  indicate  that  it  is  quite  likely  to  suffice  for  most  work.” 
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7.1 


PROPAGATION  OF  ERRORS  INTO  POSITION 


Since  most  of  the  measured  quantities  are  angle,  direction  cosines,  Him  coordinates,  range  sums, 
and  range  differences,  the  effects  of  errors  in  these  observations  are  present  in  rectangular  coordinates 
derived  from  these  quantities.  While  it  is  most  convenient  to  consider  the  errors  in  the  actual 
measurements  of  the  instrumentation  system  when  discussing  its  accuracy,  these  errors  do  not 
necessarily  give  the  desired  information.  The  relation  between  errors  in  range,  azimuth,  and  elevation 
and  errors  in  derived  position  data  depends  upon  three  factors:  (1)  the  location  and  geometry  of  the 
instrumentation  sites,  (2)  the  errors  in  the  measurements  of  the  system,  and  (3)  the  position  of  the 
point  in  space  that  is  to  be  determined. 

Methods  can  be  developed  and  programmed  which  give  estimates  of  the  errors  in  rectangular 
coordinates  when  items  (1),  (2),  (3)  in  the  above  paragraph  are  taken  into  consideration.  The  end 
product  of  these  programs  tests  estimates  of  errors  in  rectangular  coordinates,  and  is  referred  to  as 
geometric  dilution  of  precision  (GDOP).  For  examining  the  capabilities  of  given  instrumentation  on  a 
particular  test,  significant  estimates  of  the  rectangular  coordinate  accuracies  can  be  obtained  by  using 
the  nominal  trajectory  for  the  test,  the  actual  instrument  sites  to  be  used,  and  the  best  estimates  of  the 
errors  of  measurement  for  the  particular  systems.  Further  information  on  GDOP  may  be  obtained 
from  the  Bibliography. 

Although  GDOP  is  widely  used,  the  more  sophisticated  concept  of  an  ellipsoid  of  error  is  gaining 
acceptance  as  spatial  measurements  are  made  with  higher  orders  of  precision.  Whereas  GDOP  utilizes 
only  information  from  the  diagonal  terms  (oj^,  o^,  o^)  of  the  variance-covariance  matrix,  the 
ellipsoid  of  error  utilizes  all  of  the  information  which  can  be  extracted  from  the  entire  matrix.  If  an 
ellipsoid  of  error  is  used,  confidence  regions  can  be  set  up  about  each  spatial  position  point  and  the 
orientation  of  the  axes  of  the  ellipsoid  can  be  obtained.  This  is  discussed  in  more  detail  in  Section  8.0. 

The  basic  problem  in  error  propagation  is  to  determine  the  effect  of  the  random  errors  of 
measurement  in  the  particular  instruments  on  the  rectangular  coordinate  data.  Suppose,  for  example, 
that  the  quantities,  a,  0,  y,  are  measured  at  time,  t,  with  errors  A a,  A0,  Ay.  If  there  is  a  known 
mathematical  relation  between  X,  Y,  Z  and  a,  0,y: 


X  =  f,  (aJ3,y), 


Y  =  (2(aJ3,y), 


(7.1.0) 


Z=f3  (0,0,7), 


then  the  errors  in  X,  Y,  and  Z  produced  by  Aa,  A0,  Ay  can  be  approximated  by 

AX  "  If  Aa  +  If  AB  +  If 
AY-f  Aa  +  f 

Az "  If Aa  +  It  AB  +  37  AY’ 
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(7.1.1) 


where  the  partial  derivatives  are  evaluated  using  f, ,  f2,  f3  and  the  measured  values  of  a,  (3,  y.  This  type 
of  estimation  of  the  errors  is  also  satisfactory  when  discussing  the  effects  of  systematic  errors  if  these 
errors  are  relatively  small  with  respect  to  the  magnitudes  of  X,  Y  and  Z. 

A  discussion  for  other  systems  involving  multiple  stations  and  least  square  solutions  is  based  upon 
the  same  ideas.  However,  they  ire  complicated  by  the  fact  that  more  variables  are  involved  and  more 
complex  mathematical  manipulation  must  be  carried  out.  For  simplicity,  methods  are  illustrated  here 
for  radar. 

The  relations  between  the  measured  quantities,  azimuth  (A),  elevation  (E),  and  range  (R)  and  the 
space  position  (X,Y,Z)  in  left  handed  coordinates  are: 

X  =  R  cos  A  cos  E 

Y  =  R  sin  A  cos  E  (7.1.2) 

Z  =  R  sin  E. 

Errors  may  exist  singly  in  one  of  the  measurement  R,  A,  E  or  simultaneously  in  R,  A,  E.  Using 
the  above  relations  and  simultaneous  errors  in  R,  A,  E,  the  First  order  effects  of  these  errors  upon  X, 
Y,  Z  are  given  by: 


AX  *  AR  cos  E  cos  A  -  R  AA  cos  E  sin  A  -  R  AE  sin  E  cos  A 

AY  ■  AR  cos  E  sin  A  +  R  AA  cos  E  cos  A  -  R  AE  sin  E  sin  A 


(7.1.3) 


AZ  i  AR  sin  E  +  R  AE  cos  E 

Using  the  covariance  equation  (7.0.10),  we  estimate  the  elements  of  the  variance-covariance 
matrix  a  for  the  radar  assuming  the  errors  in  R,  A,  E  are  uncorrelated.  The  estimates  are 
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■  o2  cos2A  cos2E  +  02R2 


cos 2 A  sin2E  +  o? R2  sin2A  cos2E, 
A 


■  a2  sin2A  cos2E  +  0^R2  cos2A  cos2E  +  02R2  sin2A  sin2E, 


o 


2 

Z 


0 


2 

E 


■  02  sin2E  +  02R2  cos2E, 


_  SIX  9Y  _2  .  9X  9Y  *  ,  M  9Y  »  (7.1.4) 

°xy  "  9R  9R  °R  9A  9A  °A  +  9E  9E  °E 

-  (cos  A  sin  A  cos2E)  02  -  (R2  cos  E2  sin  A  cos  A)  0^ 

+  (R2  cos  A  sln2E  sin  A)  C2, 


0 


xz 


9X  3Z  „i  4.  M  JZ  „2  .  9X  9Z  2 
9R  9R  R  9A  9A  A  9E  9E  E 


*  (cos  A  cos  E  sin  E)  02 


(R2  cos  A  sin  E  cos  E)  02, 

& 


9Y  12.  -a  .  H  12,  2  .  il  9z  2 
9R  9R  aR  9A  9A  A  9E  9E  E 


■  (sin  A  cos  E  sin  E)  02 


(R2  sin  A  sin  E  cos  E)  02. 

£ 
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The  variances  and  covariances  correspond  to  the  error  distribution  which  has  been  propagated 
into  a  left  handed  coordinate  system  tangent  to  the  spheroid  with  the  radar  at  the  origin  of  the 
system. 

In  many  instances,  it  is  desired  to  propagate  the  error  into  a  coordinate  system  with  a  different 
origin  and  orientation.  For  example,  it  is  often  desired  to  propagate  the  error  into  a  coordinate  system 
whose  origin  is  a  launching  pad  and  whose  +X  axis  is  down  the  launch  azimuth.  In  this  case  the  R  A  E 
data  is  first  transformed  into  a  local  tangent  plane  system,  (Xjp.Y-ppZjp).  Based  on  the  geodetic 
position  (0,'X)  of  the  radar,  the  tangent  plane  coordinates  are  transformed  so  that  each  axis  is  parallel 
to  a  geocentric  (earth  centered)  system,  translated  to  geocentric  origin,  translated  to  new  origin 
(launch  pad),  then  the  axes  are  rotated  into  the  tangent  plane  system  with  the  +X  axis  parallel  to  the 
desired  orientation.  This  process  is  accomplished  by  the  following  matrix  equation 


*  “V-'V 


where 

[R]  *  matrix  which  rotates  radar  tangent  plane  into  geocentric  system 
[GgJ  and  [Gp]  are  geocentric  coordinates  of  radar  and  launch  pad 

[P]  ^  =  Matrix  which  rotates  axes  parallel  to  local  tangent  plane  so  translation  from  earth’s  center 
can  be  obtained.  [P]  corresponds  to  [R]  only  is  for  the  launch  pad. 

a  =  orientation  azimuth  of  transformed  system 


LZTPJa 


Yjp  =  local  tangent  plane  coordinates  with  origin  at  launch  pad  P  with  +X  axis 
oriented  parallel  to  a  degrees 


The  matrices  [R]  and  [P]  are  linear  transformations  whose  elements  are  constants  which  are 
functions  of  a,  <p  and  X.  Since  the  transformations  are  linear,  the  propagated  estimate  for  the 
covariance  matrix  in  the  new  coordinate  system  is 


«p  - A  V  * 


(7.1.6) 


where  ■  covariance  matrix  with  respect  to  radar  tangent  plane 


A  -  [PK  [R] 
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It  is  pointed  out  that  the  error  propagation  techniques  can  propagate  from  X,Y,Z  back  to  the 
original  observations  of  R,A,E.  In  order  to  do  this  the  geometry  (X,Y,Z)  and  the  covariance  matrix 
must  be  known.  For  example,  for 


R  =  f(X,Y  ,Z) 


then  "  (3^)  ax  +  W  aY  +  W  °Z 


r3R 


3R> 


(7.1.7) 


r3R^  r3R 


3R^  r3R^ 


Should  it  be  necessary  to  differentiate  the  position  to  acquire  velocity,  it  is  assumed  that  any 
low-frequency  systematic  error  is  virtually  eliminated  in  the  process.  It  is  noted,  however,  that  cyclic 
errors  may  not  be  disposed  of  in  this  manner. 

It  is  noted  that,  under  certain  conditions,  the  effects  of  total  errors  can  be  propagated  into  spatial 
position  error.  If  the  variance  of  the  systematic  error  corresponding  to  the  same  point  in  space  is 
estimated  over  several  tests,  then  the  effects  of  this  can  be  propagated  into  spatial  position  error.  In 
this  case  the  variance  of  the  total  error  input  for  propagation  would  be 

2  2  2 

0 TOTAL  "  ^SYSTEMATIC  +  0 RANDOM  (7.1.8) 


and  the  propagated  error  would  have  a  variance-covariance  matrix  corresponding  to  total  spatial 
position  error. 


7.2  PROPAGATION  OF  ERRORS  INTO  VELOCITY  AND  ACCELERATION 


In  general,  position  and  range  rate  data  are  the  basic  measurements  received  by  instrumentation 
systems  which  are  external  to  a  vehicle  in  flight.  The  basic  trajectory-related  measurements  from  on 
board  systems  are  acceleration  components  parallel  to  the  missile  coordinate  system.  In  the  first  case, 
velocity  and  acceleration  data  are  obtained  by  numerically  estimating  the  first  and  second  derivatives 
of  the  position  data.  In  the  latter,  velocity  and  position  are  obtained  by  successive  integration  with 
respect  to  time. 

Very  low  frequency  systematic  errors  in  position  data  may  be  regarded  as  constant  over  short 
spans  of  time.  In  this  instance  a  measured  variable  can  be  represented  as  a  function  of  time  by 

V  =  F(t)  +  B  (7.2.1) 

where  B  is  the  bias  error.  Upon  numerical  differentiation  with  respect  to  time  we  obtain 

V  =  F'(t)  (7.2.2) 

which  is  independent  of  B.  The  significance  of  this  fact  is  that  a  differentiation  derived  from  a 
measured  function  with  low  frequency  bias  error  is  virtually  free  of  this  bias.  However,  higher 
frequency  cyclic  error  and  random  error  will  persist  and  must  be  considered.  The  amount  of  random 
error  for  a  particular  instant  of  time  is,  of  course,  unknown  and  cannot  be  predicted  as  a  function  of 
time,  but  under  the  assumption  that  it  is  a  stationary  process  (at  least  during  short  spans  of  data),  its 
statistical  properties  remain  constant. 

Inertial  guidance  systems  currently  in  use  in  many  missile  programs  provide  examples  of 
acceleration  measuring  devices.  The  functions  these  instruments  perform  are  varied  and  complicated. 
Inertial  guidance  systems  are  discussed  in  more  detail  in  Section  10.0.  In  this  section  we  are  concerned 
primarily  with  errors  propagated  due  to  numerical  differentiation  of  position  data  with  respect  to 
time. 


If  “raw”  velocities  in  space  position  data  were  computed  by  successive  differences,  the  “average” 
velocity  estimate  over  the  time  increment  would  be 

X  ;  AX2i+l  "  At  (Xi  “  Xi-iK  (7-2.3) 

2 

The  Xj’s  are  random  variables  and  assumed  to  be  independent  with  mean  n  and  variance  a^. 
Then,  using  the  propagation  techniques  of  Section  7.0,  the  variance  in  AX  is 


Var 


+  Ot 


X± 11 


At' 


(7.2.4) 


Preceding  page  blank 
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From  (7.2.4)  it  is  obvious  that  the  random  variable  AX  contains  more  random  error  than  X.  For 
example,  if  c^=  2  feet  and  At  =  0.1  second,  then  the  error  in  raw  velocity  would  be  ±20j2or  ±28.28 
feet  per  second.  In  the  case  of  “raw”  acceleration  it  is  much  worse,  and  for  a ^  =  2  feet,  ojj  =  ±500 
feet  per  second2 !  These  examples  were  given  to  show  that  filtering/smoothing  techniques  are  needed 
to  achieve  reasonable  X  and  X  estimates  when  derived  from  numerically  differentiating  position  data. 

Most  methods  presently  used  filter  the  position  data  as  a  function  of  time  and  then  numerically 
differentiate  the  filtered  data.  The  most  common  type  of  filter  used  is  a  kth  degree  polynomial  which 
is  fitted  by  some  criterion  (usually  least  squares)  to  a  span  of  m  points.  The  first  and  second 
derivatives  of  the  polynomial  are  evaluated  at  an  appropriate  time,  usually  the  midpoint,  since  this 
gives  velocity  and  acceleration  data  with  smaller  random  errors  than  end  point  estimates.  The  point  of 
evaluation  which  gives  minimum  random  error  is  between  the  end  and  center  points.  The  moving  arc 
technique  of  fitting  the  polynomial  of  kth  degree  of  the  data,  X2  to  Xm+}  and  t2  to  tm+i ,  and 
evaluating  first  and  second  derivatives  is  used.  This  is  continued  until  all  desired  velocity  and 
acceleration  components  are  computed.  In  the  so-called  simple  differencing  method,  the  degree  of  the 
polynomial  is  2.  If  the  data  is  equally  spaced  in  time,  i.e., 


tj  =  t0  +  (At)  j. 


(7.2.5) 


then  the  polynomial  fdters  for  velocity  and  acceleration  will  consist  of  weights  bj  and  Cj  and  the 
velocity  and  acceleration  estimates  are  made  from  the  following  linear  combinations: 


m 

At  ^  bj  Xi-j+d 
J-l 


and 


(7.2.6) 


m 


J-l 


'J  Ai-j+d 


(7.2.7) 


respectively,  where  d  is  related  to  the  delay  (i.e.,  if  m  =  2n+l  and  evaluation  is  at  the  end  point  then  d 
=  m).  The  bj  and  cj  depend  on  the  degree  k  of  the  polynomial,  the  span  of  points,  m,  used  in  the 
moving  average,  ana  the  number  d.  The  bj  and  Cj  are  fixed  constants  and  are  precomputed.  Using 
formulas  of  the  type  above  and  assuming  that  a  polynomial  of  kth  degree  approximates  the  data  to  a 
sufficient  degree  of  accuracy  over  the  span  of  m  points,  and  that  the  errors  in  successive  values  of  Xj 
are  uncorrelated,  the  relation  between  the  random  errors  in  position  data  and  the  random  errors  in 
velocity  and  acceleration  data  can  be  expressed  in  the  form 

q  m 

»j)  (7-2-8> 

J-l 
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and 


J-l 


Similar  relations  can  be  obtained  if  the  errors  are  correlated  and  the  form  of  the  correlation  is  known. 
In  the  case  of  midpoint  smoothing  the  equations  for  and  o^become 

A  8X  /  12  =  M  (7.2.10) 

°Xm  It  J  o(nz  -  1)  Atm  * 


for  second  degree  polynomial,  and 


a 


x 


°X_  /  720 

(At)2  v  m(mz  -  1)  (mz  -  4) 


/720 

r* 

✓in 


(7.2.11) 


for  a  third  degree  polynomial.  It  is  evident  from  the  approximations  that  if  the  sampling  rate  is 
stepped  up  and  the  smoothing  time  span  I  Atm)  held  constant,  the  velocity  and  acceleration  errors  are 
reduced  approximately  by  a  factor  of  m'/2.  Truncation  errors  limit  the  extent  to  which  this  can  be 
applied.  If  smoothing  is  performed  at  some  point  other  than  the  midpoint  the  same  general  argument 

applies. 

In  the  discussion  of  the  effects  of  random  errors  on  velocity  and  acceleration,  it  was  pointed  out 
that  the  method  presented  was  valid  only  if  the  errors  in  successive  values  of  a  given  coordinate  are 
uncorrelated.  This  assumption  is  not  true  for  a  systematic  error  by  definition.  Hence,  different 
techniques  are  required  for  determination  of  the  effects  of  systematic  errors. 

In  the  case  of  single  station  systems,  the  errors  in  velocity  and  acceleration  in  Cartesian 
coordinates  may  be  found  by  differentiating  the  XYZ  equations  with  respect  to  time  and  utilizing  the 
Gaussian  approximation. 


(7.2.12) 
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►<  *K3  X 


where  fj  =  X,  Y,  Z,  X,  Y,  2  and  CKj  =  R,  A,  E,  R,  A,  t  for  radar.  For  example,  the  variances  of  the 
velocity  errors  for  radar  become 


■  (-E  sin  E  cos  A  -  A  cos  E  sin  A) 2  0^ 

+  (-R  sin  E  cos  A  -  E  R  cos  E  cos  A  +  A  R  sin  E  sin  A)  2  02  (7.2.13) 

+  (-R  cos  E  sin  A  +  E  R  sin  E  sin  A  -  A  R  cos  A  cos  E) 2 

+  (cos  A  cos  E)2  +  (-R  cos  A  sin  E) 2  o|  +  (-R  sin  A  cos  E) 2 

*  (-E  sin  E  sin  A  +  A  cos  E  cos  A) 2  a2 

K 

+  (R  cos  E  cos  A  -  E  R  sin  E  cos  A  -  A  R  cos  E  sin  A)  2  02  (7.2.14) 

A 

•  # 

+  (-R  sin  E  sin  A  -  E  R  cos  E  sin  A  -  A  R  sin  E  cos  A)2  02 

E 

+  (sin  A  cos  E) 2  0-  +  (-R  sin  A  sin  E) 2  0*  +  (R  cos  E  cos  A) 2  0* 

R  E  A 


•  (E  cos  E)2  02  +  (R  cos  E  -  E  R  sin  E)2  a2 
Z  R  E 

+  (sin  E)2  a^+  (R  cos  E)2  a| 


(7.2.15) 
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The  expressions  for  the  covariance  terms  o^y,  O”^,  may  be  written  according  to  the 
following  expression  for  O” y. 


9Y  _2  ,  9X  3Y  2  .  3X  3Y  2 

3R  3R  °R  3A  3A  A  3E  3E  E 


3X  3Y 
3R  3R 


°R 


+ 


3X  3Y  ,  . 
3A  ft  A 


3X  3Y  2 
3i  3i  °i 


(7.2.16) 


In  a  similar  manner,  a  partitioned  variance-covariance  matrix  can  be  estimated  giving  all  position, 
velocity,  acceleration  error  variances  and  covariances.  The  matrix  would  be 
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In  the  case  of  multiple  station  systems  requiring  least  squares  solutions  for  data  reduction,  the 
variances  cannot  be  readily  determined  analytically  as  for  radar.  The  same  reasonable  assumption  is 
made,  however,  that  the  variances  of  the  instrument  errors  represents  a  normal  distribution  at  the 
same  point  in  the  trajectory  over  many  tests.  With  this  the  Gaussian  approximation  is  again  utilized, 
except  that  the  partial  derivatives  are  replaced  by  increments,  i.e., 


The  ratio  of  Afj/AOj  is  determined  by  numerically  computing  Afj  with  a  predetermined  value  of  Aotj 
with  all  other  j  values  of  the  instrument  variables  held  constant. 
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8.0  ERROR  TOLERANCE  REGIONS 


Confidence  and  tolerance  intervals  on  the  point  estimates  of  error  parameters  were  discussed  in 
Section  5.2.  The  next  points  in  question  are  confidence  regions  on  two-  and  three-dimensional  spatial 
position  data.  In  general,  the  type  of  tolerance  regions  used  consist  of  ellipses  in  the  two-dimensional 
case  and  ellipsoids  in  the  three-dimensional  cases.  (Special  cases  of  these  are  circles  and  spheres.) 


In  obtaining  any  type  of  confidence/tolerance  region  one  must  first  make  an  assumption 
concerning  the  underlying  probability  distribution  function,  and  then  know  or  have  a  good  estimate 
on  the  parameters  which  characterize  the  distribution.  The  usual  assumption  is  that  the  underlying 
distribution  is  approximately  Gaussian  which  means  that  all  nth  dimension  distributions  can  be 
characterized  in  terms  of  the  first  and  second  moments  of  the  distribution  (i.e.,  the  means  and 
variances).  Thus,  the  basic  tool  for  the  confidence  region  is  the  variance  or  variance-covariance  matrix. 
The  probability  density  function  for  the  Gaussian  distribution  in  three  dimensions  is 


f(X,Y,Z) 


7 77717771 

(2tr)  0 


-1/2  Q"1(X,YlZ) 
e 


where 


|o|  =  determinant  of  covariance  matrix  a 

Q"*  (X,Y,Z)  =  quadratic  form  of  the  inverse  covariance  matrix  o'1 
The  quadratic  form  of  the  matrix  is 


(X,  Y,  Z)  [o"1] 


-  X2o2  +  Y2o2  +  Z2o2'  +  2XYo  '  + 
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xz 


+  2YZo 

yz 


(8.0.1) 


where  the  prime  designates  the  elements  of  a"1.  The  above  relation  is  the  general  equation  of  an 
ellipsoid  and  it  can  be  transformed  into  standard  form 


1 


(8.0.2) 


by  rotating  the  axes  of  the  general  form  equation  parallel  to  the  (X,Y,Z)  coordinate  system.  It  can  be 
shown  that  this  process  is  equivalent  to  diagonalizing  the  inverse  covariance  matrix  which,  in  turn,  is 
equivalent  to  finding  the  eigenvalues  of  o'* . 

The  diagonal  matrix  D  of  o"1  is  the  matrix  of  eigenvalues 
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Preceding  pege  blank 
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(8.0.3) 


It  is  noted  that  if  the  errors  in  X,  Y  and  Z  are  independent,  then  the  general  equation  of  the  ellipsoid 
has  no  cross  product  terms  and  corresponds  to  a  matrix  a"1 ,  which  is  already  a  diagonal  matrix.  If  the 
errors  are  independent  and  a*  =  °y  =  a\  »  the  distribution  is  the  spherical  normal  distribution  (3 
dimensions)  or  the  circular  normal  when  two  axes  are  equal  and  the  third  is  zero. 

When  the  errors  in  X,  Y  and  Z  are  correlated,  then  a  has  off-diagonal  terms  and  the 
corresponding  error  ellipsoid  is  not  in  standard  form,  but  is  skewed  in  space.  The  eigenvalues  of  o'1 
provide  the  length  of  the  axes  and  the  eigenvectors  of  a"1  will  tell  the  orientation  of  each  of  the  axes 
with  respect  to  the  coordinate  system  (X,Y,Z). 
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8.1  TOLERANCE  ELLIPSOIDS 


Once  an  error  ellipsoid  has  been  obtained,  the  question  of  how  much  of  the  spatial  error 
distribution  is  contained  in  the  ellipsoid  is  encountered.  Since  the  errors  are  assumed  to  be  normal,  it 
can  be  shown  that  the  axes  of  the  standard  ellipsoid  equation  have  the  Chi-square  distribution  with 
three  degrees  of  freedom.  With  three  degrees  of  freedom,  the  probability,  p,  that  a  point  will  be  in  an 
ellipsoid  with  constant  term  C2  is 


P  - 


(8.1.0) 


The  value  of  C2  corresponding  to  p  can  be  obtained  from  any  Chi-square  table.  The  following 
table  gives  values  of  C  versus  p. 


p 

C 

.95 

2.79 

.90 

2.50 

.80 

2.15 

.70 

1.91 

.50 

1.54 

.20 

1.00 

TABLE  1 


Thus,  for  a  90%  error  tolerance  ellipse,  C  -  2.50  so  that  the  axes  would  be 

2.5 


2.5 

/M 


2.5 

*^3 


(8.1.1) 


The  values  of  C  =  1.538  and  C  =  4.0  are  the  values  of  C  corresponding  to  “probable”  ellipsoid  of 
error  and  ellipsoidal  near-certainty  error,  respectively. 

In  the  two-dimensional  case,  the  same  procedure  is  used  except  that  o'*  is  a  2x2  matrix  and  the 
standard  equation  of  the  ellipse  is  Chi-square  with  two  degrees  of  freedom.  A  table  of  p  versus  C  is 
given  below. 


P 

C 

.99 

3.37 

0.95 

2.45 

0.90 

2.146 

0.865 

2.00 

0.50 

1.177 

0.394 

1.00 

TABLE  2 
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When  C  =  1.774  we  have  a  tolerance  ellipse  called  the  probable  ellipse  of  error.  When  the  axes  are 
equal,  the  tolerance  region  corresponding  to  C  =  1.774  is  called  “Circular  Probable  Error”  (CPE). 


In  actual  practice  a*  ^  °y  ^  °z  anc^  errors  not  be  independent,  so  it  is  usually  more 
valid  to  work  with  ellipsoids  and  ellipses  rather  than  circles  or  spheres.  There  are  certain  instances, 
however,  when  it  is  desired  to  obtain  a  circular  or  spherical  tolerance  region  from  an  ellipse  or 
ellipsoid.  That  is,  it  is  desired  to  find  the  radius  of  a  circle  (or  sphere)  which  contains  the  same  number 
of  points  as  the  tolerance  ellipse  (ellipsoid)  with  axes  C/J  Xt ,  CjJ\2-  To  accomplish  this,  we  make  the 
change  of  variables  letting  x  =  r  cos  6,  y  =  r  sin  6  then  integrating  the  bivariate  elliptical  normal 
distribution  with  respect  to  r  and  9  we  find  p  as  a  function  of 


A 


Max 


and  B 


R 

__EL 


where  R  is  the  radius  corresponding  to  p.  The  following  table  gives  A  and  B  versus  p  =  0.90. 


A 

B 

1.0 

2.1460 

1.5 

2.2501 

2.0 

2.4565 

•  2.5 

2.6865 

3.0 

2.9118 

4.0 

3.3290 

5.0 

3.7058 

10.0 

5.2111 

15.0 

6.3756 

20.0 

7.3594 

50.0 

11.631 

100.0  I  16.449 


TABLE  3 


For  example,  if  the  maximum  ratio  of  o^to  a  -  is  A  =  2.0  then  the  radius  of  a  circle  containing  90%  of 
the  error  distribution  would  be  R  gg  =  2Ab(>ijo^o^. 

The  radius  R_  of  a  sphere  containing  p%  of  an  error  distribution  can  be  obtained  from  a 
probability  ellipsoid  in  a  similar  but  more  complicated  manner,  and  will  not  be  discussed  in  this  report 
(See  Reference  29  for  details). 
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9.0  DIGITAL  FILTERING/SMOOTHING  OF  DATA 


As  stated  in  Section  2.0,  the  measured  or  raw  data  contains  signal  with  error  superimposed  on  it. 
The  basic  reason  for  using  a  digital  filter  is  to  separate  or  suppress  the  noise  error  from  the  signal. 
There  are  many  different  filtering  techniques  used  at  the  various  ranges  but,  in  general,  a  digital  filter 
is  evaluated  by  the  following  criteria: 

1)  How  effectively  the  noise  components  of  the  input  data  are  attenuated  or  removed. 

2)  How  much  distortion  in  the  signal. 

3)  How  effectively  the  filter  recovers  and  smooths  the  first,  second,  and  higher  order 
derivatives. 

4)  Amount  of  serial  correlation  in  the  data. 

5)  Amount  of  computing  time  required. 

It  would  be  desirable  to  select  a  filter  which  would  give  the  best  of  each  item  above,  but 
unfortunately  this  is  not  usually  the  case.  Often  a  filter  may  be  designed  to  accomplish  1,  2,  or  3  in 
some  “best”  manner,  but  usually  Item  5  is  then  sacrificed. 

In  general,  most  filters  used  are  linear  operators  on  the  raw  data.  If  Xj  ...,Xn  is  a  set  of  raw  data, 
the  filtered  data  is  usually  expressed  as  a  linear  combination  of  subsets  or  data  spans  of  the  set.  For 
example,  if  the  filter  gives  a  point  estimate  corresponding  to  the  center  time  point  of  a  span  consisting 
of  2N+1  points  then  the  expression  for  the  output  of  the  filter  would  be 

i«+N 

*t-  l  Vt+l  (9.0.1) 

i— N 


where  the  Wj  are  filter  weights.  These  weights  are  constrained  so  that 

l  Kt  -  1  <9-°-2> 


This  is  necessary  since  the  signals  occuring  in  missile  trajectory  work  usually  have  a  very  large 
low-frequency  trend,  and  the  response  of  the  filter  at  these  low  frequencies  should  be  exactly  one  (i.e., 
no  distortion  or  biasing  at  zero  frequency). 

If  the  filter  is  an  end  point  estimation  or  “predictor”  filter,  then  the  weighted  sum  would  consist 
of  k  previous  points. 


Many  of  the  digital  filters  being  used  are  of  the  general-purpose  type.  The  parameters  which 
define  the  weights  of  such  a  filter  are  adjusted  according  to  the  sampling  rate  of  the  input  data  and 
some  estimate  as  to  the  frequency  composition  of  the  signal.  Some  types  of  general-purpose  filters  are 

1 )  Straight  moving  average 

2)  Classical  least  squares  polynomial  fit 

3)  Constrained  least  squares  polynomial  fit 

4)  Variable  span 

5)  Frequency  constraining 

The  straight  moving  average  is  comparable  to  a  center  point  estimate  from  a  first  degree 
polynomial  in  item  (2).  The  polynomial  curve  fitting  technique  gives  a  single  point  estimate  from  a 
polynomial  curve  fit  to  n  consecutive  data  points.  The  estimate  may  be  end  point,  center  point,  or  any 
other  point  on  the  curve.  The  constrained  least  squares  uses  the  classical  approach  but  adds  constraint 
conditions  on  the  filter  parameters  based  on  some  known  information  concerning  the  data  preceding 
the  current  filter  span.  The  variable  span  filter  is  any  of  the  filter  types  with  the  added  characteristic 
that  the  filter  span  is  adjusted  (“opened  up”  or  “closed  out”  and  “reinitialized”)  based  on  some 
characteristics  of  the  input  data. 

All  linear  digital  filters  are  frequency -constraining  in  that  they  suppress  certain  frequencies,  pass 
other  frequencies  without  distorting,  and  may  distort  or  amplify  other  frequencies. 

Because  of  the  low  frequency  trend  characteristics  in  trajectory  data,  the  filters  used  to  smooth 
this  data  are  of  the  “low  pass”  type.  That  is,  the  filter  weights  are  designed  to  suppress  high 
frequencies  and  to  output  the  low  frequency  components  with  a  minimum  amount  of  distortion.  In 
many  filters  such  as  the  least  square  polynomial  type,  the  cut-off  frequency  is  determined  as  a 
function  of  the  number  of  points  in  the  filter  span.  As  a  number  of  points  in  the  filter  span  becomes 
larger,  the  cut  off  frequency  becomes  lower,  i.e.,  more  smoothing  is  applied  to  the  data.  There  are 
other  filters,  however,  which  construct  the  weights  for  a  specific  cut-off  frequency,  and  the  increase  in 
the  filter  span  does  not  change  the  cut-off  frequency  of  the  filter,  but  rather  improves  the  filter’s 
frequency  response  corresponding  to  the  desired  cut-off  frequency. 

In  addition  to  low-pass  filters,  weights  can  be  constructed  for  high-pass,  band-pass,  and 
band-reject  filters.  All  filters  mentioned  thus  far  can  be  evaluated  by  estimating  their  frequency 
response  function  H(f).  The  frequency  response  is  estimated  by  taking  the  Fourier  Transform  of  the 
weights  Wj  as  a  function  of  i.  A  typical  frequency  response  versus  the  ideal  frequency  response  or  step 
function  is  given  on  the  following  page : 


A 
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(1.0) 


An  example  of  some  frequency  response  functions  was  discussed  in  Section  4.6  and  illustrated  in 
Figure  5.  The  frequency  response  function  is  very  important  when  it  comes  to  error  analysis,  because 
knowledge  concerning  the  cut-off  frequency,  the  high-frequency  “leakage”  (side  lobe  effect),  the 
sharpness  of  filter  roll-off  at  the  cut-off  frequency,  and  the  distortion  or  amplification  of  the  low 
frequencies,  can  indicate  whether  or  not  the  appropriate  filter  was  used  for  a  given  set  of  data.  The 
frequency  response  for  a  finite  number  of  weights  will  never  follow  the  ideal  response  curve.  The 
frequency  response  oscillates  or  has  side  lobes  at  the  higher  frequencies  because  the  filter  is  truncated 
in  the  time  domain,  and  significant  weights  Wj  have  been  discarded,  causing  discontinuities  in  the  time 
response  of  the  filter.  These  time  discontinuities  in  the  time  domain  result  in  the  oscillations  in  the 
frequency  domain  at  the  higher  frequencies.  To  avoid  this  result,  a  function  must  be  chosen  which 
decays  very  quickly  in  the  time  domain,  so  that  truncation  at  a  reasonable  span  length  will  discard 
data  multipliers  of  much  smaller  magnitudes.  Slowness  of  decay  in  the  time  domain  is  caused  by 
discontinuities  in  the  derivatives,  particularly  the  zeroth  and  other  low  order  derivatives,  in  the 
frequency  domain.  This  suggests  that,  where  possible,  functions  having  all  continuous  derivatives  be 
chosen  for  use. 

The  error  introduced  into  the  filtering  process  by  the  oscillations  in  the  transfer  function  (curve) 
is  inversely  proportional  to  the  number  of  points  used  in  the  filter  and  the  roll-off  frequency  specified. 
However,  there  are  means  whereby  this  response  can  be  controlled  and  a  fixed  cut  off  specified,  and  at 
the  same  time  have  a  rapid  roll  off  and  very  low  passage  in  high  frequency  areas.  Such  filters  are 
described  in  Reference  3. 

A  discussion  of  differentiation  and  prediction  filters  (see  Reference  4)  would  be  too  lengthy  to 
pursue  in  detail.  Nevertheless,  a  few  guiding  comments  will  be  added  here.  Two  frequency  filters  with 
different  cut-off  frequencies  can  be  “cascaded”  in  the  data  reduction  process  to  achieve  the  digital 
equivalent  of  cither  a  band-pass  or  a  band-elimination  filter,  depending  upon  the  arrangement  of  the 
data  multipliers  and  of  the  filter  outputs  and  residuals.  To  preserve  the  high  rate  of  change  in  slope 
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and  avoid  filter  lags  at  the  acceleration  discontinuities  at  staging,  it  becomes  necessary  to  have  a 
'ariable  span  filter,  and  collapse  the  filter  weighting  function  into  the  discontinuities  and  build  up  the 
weighting  function  again  on  the  other  side.  Unfortunately,  such  a  process  does  not,  in  general,  preserve 
the  desired  qualities  of  the  function.  Therefore,  if  an  integration  step  accompanies  the  smoothing 
process,  such  as  occurs  in  deriving  velocity  from  acceleration,  it  is  necessary  either  to  filter  through  the 
discontinuity,  thereby  sacrificing  the  high-frequency  response  of  the  data,  or  to  perform  the 
integration  process  first  and  then  filter,  collapsing  the  filter  when  desired.  However,  this  does  require 
additional  computations  since  the  raw  data  is  generally  at  a  higher  sample  rate  than  the  smoothed 
data,  and  in  common  with  the  prefiltering  approach,  careful  editing  of  the  raw  data  is  required  to 
prevent  integration  errors  due  to  large  noise  spikes  in  the  raw  data  (see  Reference  2). 

Another  type  of  frequency  digital  filtering  has  become  feasible  and  desirable  due  to  a  recently 
derived  algorithm  for  spectral  analysis  called  the  Fast  Fourier  Transform  (FFT)  (Reference  27). 
Besides  making  digital  spectrum  analysis  more  attractive  from  an  economic  standpoint,  the  FFT  has 
enabled  many  to  change  their  concepts  of  digital  filtering,  in  that  the  intellectually  appealing  approach 
of  filtering  in  the  frequency  domain  is  now  often  simpler,  faster,  and  just  as  effective  as  filtering  in  the 
time  domain,  although  two  transforms  between  the  time  and  frequency  domain  are  employed  in  the 
process.  Suppose,  for  example,  we  desire  to  filter  the  data  Xj  ,X2,...,Xn  with  a  low-pass  filter  at  some 
cut-off  frequency  fc<  Instead  of  computing  the  weights  for  a  time  domain  convolution  we  use  the  FFT 
to  estimate  the  real  and  imaginary  parts  of  the  spectrum  of  the  time  series.  Once  we  obtain  the 
spectrum  we  truncate  it  at  the  appropriate  cut-off  frequency  then  take  the  inverse  FFT.  The  resulting 
data  in  the  time  domain  corresponds  to  the  truncated  spectrum  and,  as  a  result,  has  been  “low  pass” 
filtered  with  cut-off  fc, 

In  addition  to  the  so-called  general-purpose  filters,  there  are  special-purpose  digital  filters  which 
obtain  each  filtered  data  point  from  the  preceding  data  point  using  the  recursive  relationships  and 
information  concerning  the  variance  of  preceding  data  contained  in  the  “state  vector.”  The  state 
vector  describes  the  mathematical  model  of  the  total  system  generating  the  data  which  is  to  be 
filtered.  For  this  filter  to  function  in  an  optimum  manner,  the  following  must  be  given:  (1)  system 
parameters  and  their  configuration,  (2)  state  vector  initial  values  and  variance-covariance  matrix  of 
state  variables,  (3)  time  and  effects  of  special  events.  The  special-purpose  filter  is  much  more 
sophisticated  from  a  mathematical  standpoint  than  the  general-purpose  or  simple  digital  filter.  The 
special-purpose  filter  can  be  a  very  effective  analysis  tool  especially  when  used  in  conjunction  with 
multi-instrument  solutions.  One  of  the  more  well  known  special-purpose  filters  is  the  Kalman  filter 
(see  Reference  28). 
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9.1 


EFFECTS  OF  SMOOTHING  ON  DATA 


There  is  a  diversity  of  opinion  among  people  who  process  and  analyze  data  concerning  the  many 
numerical  methods  used  in  smoothing  measured  quantities  and  in  the  retrieval  of  information  from 
erroneous  observations.  The  literature  on  this  subject  is  voluminous.  Of  the  many  methods,  all  use 
assumptions  on  the  functional  form  of  the  basic  data  trend  or  the  statistical  properties  or  origin  of 
errors  in  an  effort  to  obtain  a  numerical  process  which  will,  within  the  imposed  constraints,  improve 
the  data  by  minimizing  errors.  The  choice  of  a  technique  must  depend  upon  the  objectives  sought. 
Some  investigation  of  a  given  technique’s  potential  in  describing  the  data  must  be  made.  For  example, 
statistical  smoothing  methods  depend  upon  the  character  of  errors,  their  distribution,  their  variances, 
and  their  dependence  or  independence.  One  might  choose  a  second-degree  moving  arc  polynomial  over 
successive  51  point  spans  of  10  samples  per  second  on  position  data  to  achieve  a  certain  reduction  in 
error  variance. 

Least  squares  polynomial  smoothing  has  long  been  used  as  a  general-purpose  filter.  How  well  the 
objectives  may  be  achieved  in  choosing  a  smoothing  technique  may  be  determined  by  the  knowledge 
of  the  functional  equations  of  motion  and  the  error  characteristics.  Figures  10  thru  17  contain 
important  facts  relative  to  polynomial  smoothing  and  Figure  18  corresponds  to  a  particular  set  of 
frequency-constraining  digital  filters.  A  description  of  each  follows: 

Figure  10  shows  estimates  of  position  and  velocity  error  using  a  polynomial  of  degree  two  and 
various  (Ns)  smoothing  intervals.  Evaluation  is  at  midpoint  of  smoothing  interval. 

Figure  1 1  gives  the  relative  smoothing  in  velocity  determinations  obtained  by  using  a  span  of 
2n+l  points  and  evaluating  at  points  other  than  the  midpoint  for  the  case  where  X  can  be  expressed  as 
a  quadratic  function  of  time  and  the  errors  in  X  are  uncorrelated. 

Figure  12  give  curves  for  higher  degree  polynomial  fits  for  both  velocity  and  accleration  errors. 

Figure  17  gives  the  values 
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for  center  point  smoothing  as  a  function  of  the  number  N  of  points  in  the  smoothing  span. 

Figure  18  gives  the  reduction  factor  in  the  standard  deviation  of  the  random  error  of  position 
data.  The  family  of  curves  represents  low-pass  filters  with  spans  from  3  to  400  points  with  cut-off 
frequencies  from  0.0200  hertz  to  0.500  hertz. 


83 


Estimates  of  Errors  in  Velocity  Data  (oy)  vs.  Errors 
Position  Data  (ap)  at  Midpoints  of  Dif.  Quadratic  Smoothing 

Figure  10 


in 

Intervals (Ng) 


Estimated  Position 


Error  op  (feet) 
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NO.  OF  POINTS  FROM  MIDPOINT 
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FIGURE  18 


VARIANCE  REDUCTION  CURVE 

FOR 

WEIGHTING  FACTOR  &  NUMBER 

FILTER  TYPE  «  LOW 

CUT-OFF  FREQ  *  .0200  .1400  .2600  .3600  .5000 


10.0 


INERTIAL  GUIDANCE  DATA 


Traditionally,  attention  and  effort  have  been  centered  upon  attaining  velocity  and  accleration 
data  from  position  data.  Techniques  for  determining  the  accuracies  of  the  derived  quantities  have  been 
discussed.  Recently,  there  has  been  a  development  of  strong  interest  in  utilizing  data  from 
accelerometers  of  guidance  systems  to  assess  velocity  accuracies.  Historically,  this  may  appear  to  be  a 
paradox.  Development  of  the  many  tracking  systems  on  the  various  ranges  had,  as  a  primary  objective, 
the  evaluation  of  guidance  systems.  Hence,  it  seems  odd  now  to  consider  the  use  of  guidance  systems 
to  evaluate  tracking  systems.  However,  the  weight  of  accumulated  evidence  and  experience  has  forced 
the  conclusion  that  the  quality,  performance,  and  reliability  of  these  systems  are  high.  Nevertheless, 
the  simple  integration  of  acceleration  data  to  produce  velocities  has  an  inherent  disadvantage  in  that 
additional  unknowns  are  introduced. 

In  view  of  the  above,  techniques  have  been  developed  for  combining  this  data  with  that  of  other 
sources  to  establish  trajectories.  It  has  been  suggested  that  errors  in  velocity  arise  from  accelerometer 
errors  such  as  bias,  scale  factor,  quadratic  nonlinearity,  platform  misalignment,  gyro  mass  unbalance, 
and  gyro  constant  drift.  Errors  from  these  sources  are  relatively  small  during  the  initial  portions  of 
flight,  but  after  50  seconds  or  so  tend  to  exceed  those  present  in  data  from  the  better  CW  Doppler 
interferometers.  For  a  typical  system,  these  errors  have  been  identified  as  follows,  in  the  order  listed: 


AK„  "  M 

K  -  30  x  10_6g 

0 

AXj  -  I^X  , 

Kj  ■  10  ppm 

AX2  -  K2/a2dt  , 

K2  -  3  ppm/g 

ax3  -  k3z 

X3  ■  2  sec. 

ax4  -  k4z2 

K4  ■  0.01  deg/hr/g 

AX5  ■  Ksa  tdt  , 

Z 

K5  ■  0.05  deg/hr. 

It  may  be  observed  that  few  of  these  are  identifiable  in  functional  form  and  in  certain  cases  a  high 
degree  of  correlation  is  present.  Furthermore,  some  of  these  are  time  dependent  and  tend  to  grow  with 
time.  A  factor  not  appearing  above  but  one  exerting  considerable  influence  on  the  analysis  of  velocity 
is  the  propagation  of  tracking  noise  into  the  velocity  vector.  This  also  tends  to  grow  with  time,  and 
hence,  tends  to  obliterate  the  data  towards  the  end  of  powered  flight.  It  should  be  noted  that  on 
conventional  solid-propellant  vehicles,  loss  of  guidance  data  is  usual  following  missile  burnout  because 
of  the  dynamics  of  the  vehicle  attendant  to  thrust  termination.  Consequently,  no  free-fall  flight  data 
are  ordinarily  available. 

For  powered  flight  data  analysis,  velocity  error  estimates  may  be  obtained  directly  by  comparing 
the  velocity  data  derived  from  the  tracking  systems  positional  data  with  inertial  guidance  velocity 
data.  Inertial  guidance  data  are  subject  to  long  term  drift  errors,  but  exhibit  good  characteristics  in  the 
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very  low  frequency  errors  (errors  below,  .05  cycles  per  second);  whereas,  the  data  from  external 
systems  are  subject  to  errors  in  the  high  frequency  and  middle  frequency  regions  (errors  above  .05 
cycles  per  second),  but  exhibit  good  characteristics  in  the  low  frequency  or  near  zero  frequency 
region.  It  is  this  very  fact,  i.e.,  that  these  differences  exist  in  the  frequency  characteristics  of  the  two 
systems,  that  the  guidance  data  may  be  used  as  a  standard  for  noise  estimation  in  the  external  data, 
and  the  external  data,  in  turn,  may  be  used  as  a  standard  for  the  estimation  of  long  term  drifts  or 
biases  in  the  guidance  data. 

The  analysis  of  free-fall  data  may  proceed  along  two  lines.  It  may  be  assumed  that  the  observed 
X,  Y,  Z  data  are  serially  correlated  and  processed  through  smoothing  filters  in  order  to  obtain  velocity 
data  (X,  V,  Z)  and  acceleration  data  (X,  V,  2).  Then  estimates  of  the  error  in  the  derivative  data  may 
be  obtained  as  a  function  of  the  variance  estimate  of  the  random  +  cyclic  error  in  position.  The 
position  error  variance  which  propagates  into  velocity  and  acceleration  is  a  composite  of  the  high  and 
low  frequency  noise  in  the  data.  This  error  is  also  referred  to  as  the  random  +  correlated  error,  or 
perhaps  most  appropriately  described  as  the  total  error  minus  the  constant  bias  error.  Constant  bias 
errors  in  the  position  data  introduce  negligible  effects  in  the  derivative  velocity  and  acceleration  data, 
since  by  their  very  nature  they  tend  to  disappear  in  the  process  of  differentiation.  It  is  (cyclic)  low 
frequency  error  which  has  the  most  deleterious  effect  on  the  quality  of  the  derivative  data.  Smoothing 
over  spans  appreciably  shorter  than  the  typical  low  frequency  does  little  to  reduce  its  influence.  The 
rigorous  propagation  of  correlated  errors  through  any  specified  smoothing  filter  requires  a  knowledge 
of  the  autocorrelation  function  of  the  errors.  If  the  autocorrelation  coefficients  for  successive  values 
of  X  and  1,  Pj ,  P2 ,  P3...,  the  variance  of  the  value  of  X  obtained  from  the  smoothing  filter 


m 


bj  xi+j 


(10.0.2) 


is  given  by 

°2X  "  °x  b  R  bT  (10.0.3) 

in  which  b  is  the  row  vector  of  data  multipliers 
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-m+1 


•  •  • 


b-lb0  V 


b 

m-l 


b  ) 

m 


(10.0.4) 
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R  is  Che  auto* 

-correlation 

matrix: 

r- 

1 

Pi 

P2 

?3 

•  •  • 

P2m 

Pi 

1 

Pi 

P2 

•  •  • 

p2m-l 

P2 

Pi 

1 

Pi 

•  •  • 

p2m-2 

R  - 

P3 

• 

P2 

• 

Pi 

• 

1 

• 

•  •  • 

p2m-3 

• 

• 

• 

p2m 

• 

• 

p2m-l 

• 

• 

p2m-2 

• 

• 

P2m-3 

• 

• 

1 

and  Is  the  variance  of 


(10.0.5) 


Similarly,  the  variance  of  X  obtained  from  the  smoothing  filter 


m 


Xi  ^  cj  x1+j 
j«-m 

is  given  by 


(10.0.6) 


9  9  X 

■  a*  c  R  c 
^  x 

in  which  c  is  the  row  vector  of  data  multipliers 

c  *  (c  c  , ,  . . .  c  ,  c  ) , 

-m  -m+i  m-1  m 


(10.0.7) 


(10.0.8) 


R  is  the  autocorrelation  matrix. 

An  estimate  of  the  autocorrelation  function  is  best  obtained  from  the  residuals  about  a  free  fall 
ellipse  (Keplerian  ellipse)  fitted  to  a  long  span  (preferably  100  to  200  seconds)  of  post  burnout  data. 
The  autocorrelation  function  may  be  obtained  from  the  power  spectrum  analysis  computer  program. 
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There  have  been  many  difficulties  associated  with  the  use  of  aircraft  and  ballistic  camera  coverage 
on  missile  tests  to  calibrate  external  ti  eking  system  velocity  errors  to  acceptable  accuracies,  which 
have  led  to  the  use  of  a  second  techniqu  .  This  method  consists  of  fitting  the  Keplerian  equations  of 
motion  to  data  taken  at  both  ends  of  the  tree  flight  trajectory,  i.e.,  data  taken  by  an  up-range  tracking 
system  of  the  missile  position  after  missile  burnout  and  data  taken  by  a  down-range  tracking  system  of 
missile  position  prior  to  atmospheric  re-entry.  Position  parameters  are  fitted  to  a  Taylor  series 
expansion  of  the  equations  of  motion  truncated  at  the  fourth  power.  Because  the  time  between  the 
two  data  spans  is  accurate,  velocities  obtained  from  the  trajectory  reconstructed  in  this  manner  will 
also  be  accurate  (typically  on  the  order  of  0.1  fps)  although  there  may  be  sizeable  position 
uncertainties  (typically  200-300  feet)  in  the  data  taken  at  each  end. 


At  t  -  tQ 


R  +  V  (T  -  t  )  + 

0  0  o 


T 

t 

0 


t 

g  dtdt. 
t 

o 


(10.0.9) 


Hence , 


t 

g  dtdt. 
t 

o  o 


(10.0.10) 
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The  derivation  of  this  result  is  not  difficult.  From  the  Figure,  R0  is  the  radius  at  injection  time  t0 
and  T  is  time  of  re-entry  where  radius  is  Kj.  Since  this  portion  of  the  trajectory  is  substantially 
outside  the  atmosphere  the  use  of  no-drag  formulation  with  only  gravity  force  acting  is  justified.  Thus, 
acceleration  is: 


a=  g 


or 


dv 

dt  "  8 
dv  -  g  dt 


V 


r  T 


g  dt  +  Cj 
t 

o 


(10.0.11) 


and 


R  -  Vq(T  -  tQ)  + 


g  dtdt  +  Cz 


t  Jt 

o  o 


or 


rT 


RT  -  Ro  +  Vo(T  -  t0)  + 


g  dtdt. 


o  o 


(10.0.12) 


This  may  then  be  solved  for  Vq. 


11.0 


ERROR  ANALYSIS  IN  THE  FREQUENCY  DOMAIN 


A  very  effective  tool  to  use  in  the  evaluation  of  error  is  spectral  analysis.  It  is  assumed  that  the 
errors  e£  are  a  realization  from  a  stochastic  process  which  is  stationary  in  the  wide  sense.  That  is, 
the  process  has  a  constant  mean  and  its  autocovariance  is  independent  of  time,  but  is  dependent  on 
time  lag.  That  is 


E(et)  =  M(t)  =  CONSTANT 


and 


R£  (k)  -  E[(et  -  u)(et+k  -  u) ] .  k  -  0,1,2 . m. 


(11.0.1) 


The  assumption  of  stationarity  is  generally  true  over  subsets  of  the  data,  and  if  et  is  not  stationary  it 
can  usually  be  "prewhitened”  or  “detrended”  so  that  nonstationary  effects  can  be  removed  from  the 
data.  The  variance  of  e£  is  obtained  from  R(k)  when  k  =  0,  i.e.,  o^2  =  R(O).  The  autocorrelation 
function  of  e£  is  given  by 


P(k) 


R£kl 

R(0) ' 


It  can  be  shown  that  if  e£  is  a  stationary  process  then  the  autocovariance  has  the  following 
properties 

R(0)  >  0 
R(-k)  -  R(k) 

| R(k) |  <R(0) 

lim  R(k)  -  0  (11.0.2) 

k-KX> 


Further,  it  can  be  shown  that  the  autocorrelation  matrix 


p(p) 

p(l) 

p(2)  ... 

p(p) 

[p] 

p(-l) 

• 

P(0) 

• 

p(l)  ... 

• 

p(p-l) 

• 

• 

• 

p(-p) 

• 

• 

p(-p+l) 

« 

• 

•  e  • 

• 

• 

P(0) 

is  positive  definite 


(11.0.3) 
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The  spectral  density  function  for  continuous  et  is 


e  iwTR(T)dT 


(11.0.4) 


and  if  is  discrete  (stochastic  sequence)  then 


.  i-  Y 

L 


-lux,,/  . 
e  R(t). 


(11.0.5) 


Conversely,  if  we  know  the  spectral  density  f(oi)  then 


R(t)  -  eluTf(ui)da) 


(11.0.6) 


ela,Tf(u))du 


(11.0.7) 


for  the  continuous  and  discrete  cases  respectively.  In  digital  analysis  it  is  clear  we  are  interested  in 
using  the  discrete  case  to  approximate  the  continuous  function,  and  from  (11.0.7)  we  can  see  that  the 
quantities 


(11.0.8) 


i--R(-T),  (t-0,+1,+2 . ), 

L  TT 
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are  the  Fourier  coefficients  of  the  function  f(oo).  In  actuality  et  is  generally  continuous  but  a 
realization  et  is  usually  discrete.  In  addition  to  the  spectral  density  function,  another  important 
function  in  the  frequency  domain  is  the  spectral  distribution  function 


F(X) 


[  f(u)du 


-n 


(11.0.9) 


which  distributes  the  variance  or  power  as  a  function  of  frequency  increments.  If  the  Fourier 
Transform  of  the  autocorrelation  function  is  taken,  then  the  total  area  under  the  spectral  density 
curve  is  unity.  If  the  Fourier  Transform  of  the  autocovariance  is  taken,  then  the  total  area  under  the 
curve  is  R(o)  =  a 2  (in  this  case  the  spectral  density  is  often  called  the  power  spectrum). 


In  actual  practice  the  autocovariance  is  approximated  from  n  error  estimates  ei  ,...,en  over  a  total 
of  m  <  n  lags. 


n-p 


R 


l  (ei-  ^‘i+p  p  -0.1.2 . « 


i-i 


where 


(11.0.10) 


■  Re(pAt),  At  ■  time  increment. 

The  “raw”  spectral  density  is  estimated  by  numerically  estimating  the  Fourier  transform  of  R^P)  using 
the  trapezoidal  rule 


f(s) 


M 


24t  j  t  R(P) 
tt  L  p  e 


P-0 


cos 


s-0,1, . . .  ,M 


where 
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P 


(11.0.11) 


f  (s) 


f(“s) 
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The  raw  spectral  estimates  are  smoothed  using  a  Hamming  spectral  window.  In  the  frequency 
domain  the  Hamming  weights  give  a  three-point  moving  average  and  the  filtered  spectrum  estimates 
are 


f(0)  -  .54  f(0)  +  .46  f(l) 

f(s)  -  .23  f(s-l)  +  .54  f(s)  +  .23  f(s+l) 

A  A 

f (M)  -  ,54  f (M)  +  .46  f(M-l) 


(11.0.12) 


The  Hamming  weights  are  used  to  reduce  the  side  band  effect  caused  by  discrete  approximation  when 
frequencies  are  not  exactly  on  the  discrete  spectral  estimate  point.  The  frequency  response  of  the 
Hamming  weights  and  others  are  well  known  and  there  is  much  literature  on  the  subject.  Since 

8  TT  r\  a 

“s  '  Mi?  2"  fs' 

then  the  spectral  density  estimates  will  be  in  increments  of  Af  -  yt.v  for  ®  ^  =  2At  •  The 

frequency  is  called  the  Nyquist  or  “folding  frequency,”  f^,  and  is  the  maximum  frequency 

which  can  be  resolved  for  a  given  sampling  rate  At.  The  Nyquist  frequency  is  called  the  folding 
frequency  because  it  can  be  shown  that  if  significant  frequencies  exist  in  the  data  which  are  higher 
than  hertz,  they  will  “fold”  about  the  Nyquist  frequency  and  appear  on  the  spectrum  as  a 

lower  frequency  between  0  and  ^t"  hertz.  This  is  illustrated  in  Figure  19  below.  For  example, 
suppose  et  is  a  continuous  (analog)  function  which  contains  a  60  cycle  “hum”  component  from  a 
power  source  and  it  is  digitized  at  At  =  .05  seconds.  =  10  cycles  =  fj^  and  the  60  cycle 

component  would  falsely  appear  as  a  zero  frequency  component  or  bias  in  et.  The  folding  of  higher 
frequencies  is  called  aliasing. 


t At  •  fn  Urn  8fn  7fn 


FREQUENCY  -* 

FIGURE  19  EFFECTS  OF  FOLDING  OR  ALIASING 
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The  concept  of  aliasing  immediately  points  to  the  fact  that  if  the  data  are  sampled  with  At  such 
that  Significant  frequencies,  erroneous  frequencies  will  appear  in  the  data.  The  effects  of  these 

frequencies  on  the  error  in  position,  velocity  and  acceleration  data  can  be  highly  significant.  If  the 
data  are  such  that  a  large  At  (low  sampling  rate)  is  required,  then  one  must  be  assured  that  significant 
frequencies  do  not  exist  above  the  Nyquist  frequency.  This  can  be  assured  by  prefiltering  the  data 
using  a  higher  sampling  rate  and  a  low  pass  filter  with  a  cut-off  frequency  equal  to  the  Nyquist 
frequency  corresponding  to  the  desired  low  sampling  rate.  This  is  sometimes  necessary  when  spectral 
analysis  on  low  frequency  error  components  is  to  be  accomplished. 


Since  the  relation  between  the  autocovariance  and  the  spectral  density  function  is  known,  we  can 
utilize  these  two  functions  in  the  analysis  of  many  types  of  data  and  for  error  analysis  specifically.  For 
this,  several  examples  of  autocovariance  and  spectral  density  functions  corresponding  to  different 
error  models  are  given  (see  Figure  20). 


Suppose  e,  ,...,en  is  a  realization  from  a  random  process  which  is  independent  (uncorrelated). 
Since  the  ej  are  discrete,  the  realization  will  have  a  low  pass  frequency  band  and  since  the  realization  is 
uncorrelated,  it  is  called  white  noise.  In  this  case 


/■o2,  k  -  0 

R<k>  (11.0.13) 

U,  k  i  0 


Then,  according  to  (11.0.2)  the  corresponding  spectral  density  is 
o2 

f(co)=^“  .  Thus  stationary  sequences  of  uncorrelated  random  errors  are  characterized  by  the  fact 

that  their  spectral  densities  are  constant  over  the  interval  oj  =  0  to  co  =  rr/At.  An  example  of  the 
estimated  correlation  of  white  noise  from  radar  data  is  given  in  Figure  21. 

Suppose  now  that  the  realization  ej  ,e2,...,en  is  correlated  with  autocovariance 

R(k)  ■  a2pk,  | e [  <  1.  (11.0.14) 


m  . 

Autocorrelation  is  R(0)  =  an<^  this  ^orm  correlation  corresponds  to  a  first  order  Markov  process, 

i.e.,  et  =  p  et.j  +TJt.  The  autocorrelation  decays  exponentially  and  from  (11.0.5)  it  can  be  shown  that 


f(w) 


(11.0.15) 


which  is  an  exponential  type  spectral  density  function  (see  Figure  20). 


Figure  20  Examples  of  some  continuous  autocorrelation  functions  and  corresponding  spectral 
density  functions.  Types  A  and  C  are  the  continuous  analogs  to  a  First  and  second  order  Markov 
Process  respectively , 
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FIGURE  21  Estimate  of  autocorrelation  function  of  independent  random  error 
(white  noise). 
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Markov  processes  can  be  estimated  by  the  so-called  autoregressive  scheme.  For  example,  the 
autocorrelation  and  the  power  spectrum  can  be  estimated  for  a  particular  set  of  error  estimates.  If  the 
shape  of  the  autocorrelation  and  power  spectrum  indicate  that  the  underlying  stochastic  process  is 
Markov,  it  can  be  approximated  by  an  autoregressive  scheme.  To  accomplish  this,  the  error  estimates 
are  staggered  in  the  following  manner 


e5 

e3 

e2 

el 

e6 
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e4 

e3 
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n-3 

n 

The  First  column  is  considered  to  be  the  dependent  stochastic  variable.  The  relation  ej  =  at  ej.j  + 
0f2 ei_2  +  a3ei-3  +  a4ei-4  +  *s  l^en  estimated  by  regression  techniques.  An  analysis  of  variance  can 
indicate  which  coefficients  are  significant,  and  then  the  process  can  be  repeated  using  the  appropriate 
order  for  the  Markov  process.  An  example  of  such  a  technique  is  illustrated  in  the  autoregressive  fit  in 
Figure  23.  The  corresponding  autocorrelation  estimate  is  in  Figure  22.  An  analysis  of  variance 
corresponding  to  the  autoregressive  scheme  is  tabled  in  Figure  24.  The  analysis  of  variance  indicates 
the  process  is  of  the  first  order  because  the  regression  coefficients  a2,  a3,  a4,  are  not  significant.  Figure 
23  gives  the  graph  of  the  approximation  of  rj j.  If  ej,...,en  is  a  realization  from  the  process 
et  =  ajet_j  ■*a2et+2+,It’Tt‘s  independent  random  variable,  then  it  can  be  shown  the  autocorrelation 
function  is  an  exponentially  damped  harmonic  of  the  form. 

pk  -  [a2/2  8in  (kV  +  /8in  W  (11.0.16) 


where 
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FIGURE  22  Estimate  of  autocorrelation  function  for  first  order  MARKOV  process 
with  p  -  a\  *  0.99.  The  estimate  at  LAG  k«l  is  0.97.  Note  the  error  in  the 
estimate  for  large  lags.  RHO  is  not  supposed  to  become  negative  for  this  parti 
cular  process. 
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FIGURE  23  Estimate  of  Markov  process  using  4th  order  autoregressive  scheme. 

The  regression  coefficients  Indicate  that  the  process  Is  first  order  with  p  -  .96848. 
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The  spectral  density  corresponding  to  p^  will  be  a  band  of  frequencies  about  toQ  =W  and  the 
general  form  of  the  spectral  density  will  be 


f  (w) 


A(ui2  4-  b2) 
W2  +  2aoi2  +  b4 


where 


b  -  As*  +  W2  (11.0.17) 

a  -  k2  -  W2 

The  continuous  case  of  damped  oscillatory  covariance  function  is 

R(t)  ■  o2e  °^TUos  cj  t  (11.0.18) 

0 

Another  example  of  random  error  would  be  the  so-called  moving  average  of 
a  sequence  of  uncorrelated  random  variable.  The  stochastic  model  is 

n 

Ej  -  l  akr)j>i»  lal  <  1»  (11.0.19) 

k-1 
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A 


where  the  fj’s  are  independently  distributed  with  zero  mean.  In  this  case  it  can  be  shown  that  in  the 
limit  as  n  -*• 00  the  autocovariance  function  is  of  the  form 


R00  - 


1  - 


(11.0.20) 


which  corresponds  to  a  spectral  density  similar  to  (11.0.15)  if  a  =  p  and  oJ  =  L 

1  -  a* 

The  most  general  type  of  correlated  noise  is  one  where  the  correlation  is  varying  and  its  graph  has 
varying  amplitude  and  periods.  In  this  case  the  stochastic  error  model  is 


I  (*£  cos  oijt  +  sin  w  t) 
1-1 


(11.0.21) 


where  the  aj  and  bj  are  independent  random  variables  with  mean  zero  and  variance  o*  +  o£  .  =*  a?. 
In  this  case  the  autocovariance  function  is  1  1 


n 

R(k)  -  l  a*  cos  Wjk 
1*1 


(11.0.22) 


and  the  spectral  density  function  is  a  sequence  of  frequency  bands  made  up  of  “impulses” 
corresponding  to  a?/Aco  at  CJj.  The  analytic  form  of  f(w)  is 


u 

f(w)  —  it  £  0^  [5  (u)  -  u)1)  +  5(o)  +  coi)] 

1-1 


(11.0.23) 


where  6  is  the  Dirac  delta  function. 


(11.0.21)  is  often  the  type  of  spectral  density  encountered  in  error  analysis.  An  example  of  (11.0.23) 
is  given  in  Figure  25. 

Analysis  of  error  in  the  frequency  domain  can  be  made  more  efficient  from  a  computational 
standpoint  if  the  fast  Fourier  transform  is  used  (Reference  27).  The  autocovariance  function  can  be 
obtained  with  less  computation  time  by  taking  the  fast  Fourier  transform,  computing  the  modulus, 
smoothing,  then  taking  the  inverse  fast  Fourier  transform. 
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FREQUENCY  CPS 

FIGURE  25  Estimate  of  the  spectral  density  function  for  a  stochastric  process 
with  random  amplitudes  and  phases.  The  ordinate  estimates  the  variance  per 
frequency  increment  Af  ■  ^  • 


Once  the  frequency  content  and  autocovariance  of  the  error  et  have  been  estimated,  one  can 

filter  the  data  to  remove  both  high  frequency  random  error  and  oscillatory  systematic  errors  of 

moderately  low  frequency  (say  0.05  hertz).  In  addition,  the  autocovariance  function  can  be  used  to 

remove  biased  estimates  of  the  variance  of  the  error  by  extended  least  squares  approximation.  The 

analysis  of  data  in  the  frequency  domain  can  be  utilized  in  conjunction  with  multiple  regression  and 

analysis  of  variance.  For  example,  suppose  the  error  is  a  stationary  process  of  the  Form  (11.0.21) 

However,  the  significant  frequencies  Wj  are  not  necessarily  harmonics  and  are  not  known.  A  spectral 

analysis  would  show  significant  frequency  components  and  the  corresponding  squared  amplitudes  (see 

Figure  25).  With  k  significant  CJj’s  selected  from  the  spectrum,  (11.0.21)  could  be  approximated  using 

multiple  regression  and  solving  for  the  aj’s  and  bj’s.  An  analysis  of  variance  would  then  enable  one  to 

determine  how  good  an  approximation  would  be  made  on  the  expected  values  of  a^,  bj  and  their 

standard  deviations  a.  ,  Ol  . 

ai  Di 

As  an  example  a  regression  analysis  was  done  in  conjunction  with  error  data  corresponding  to  the 
spectrum  estimates  in  Figure  25.  The  results  of  this  analysis  are  given  in  graphical  form  in  Figure  26. 
An  analysis  of  variance  indicated  the  regression  analysis  explained  82%  of  the  total  variation  in  the 
data.  Finally,  an  estimate  of  the  spectral  distribution  function  corresponding  to  the  spectral  density 
function  is  given  in  Figure  27.  The  spectral  distribution  gives  the  estimate  of  variance  as  a  function  of 
frequency  bands.  From  the  graph  we  see  that  the  total  variance  estimate  corresponds  to  the 
intercept:  F  (0)  -  F  (^-)  =  0.954.  The  variance  estimate  over  the  band  coA  <  w3  <  ojg  is  F(cdg)  - 

F(wa)  or  =  0.0365. 
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FREQUENCY  CPS 


FIGURE  27  Estimate  of  spectral  distribution  function,  F,  corresponding  to 
spectral  density  function  in  FIGURE  25. 
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